Abstract

The drift of temperature measurements by semiconductor negative temperature coefficient thermistors is a well-known problem. This study analyzes the drift characteristics of the thermistors designed and used at the Royal Netherlands Institute for Sea Research for measuring high-frequency temperature fluctuations in the ocean. These thermistors can be calibrated to high precision and accuracy (better than 1 mK) and have very low noise levels. The thermistors can measure independently for long periods of time (more than one year), and the identification and compensation of the drift are thus essential processing steps. A laboratory analysis showing that the drift is similar, in its functional form, to the drift of commercial thermistors described in the literature is presented. An effective procedure to estimate this drift from ocean observations is described and tested using three datasets from the deep Atlantic Ocean. Since the functional form of the drift rate is, with good approximation, universal among different sensors, the procedure could easily be adapted to other datasets and, the authors argue, to measurements from thermistors by other manufacturers too.

1. Introduction

Measurements of temperature in the ocean have a centuries-old tradition. Today, temperature measurements remain important for monitoring the mean state of the ocean and its variability at different scales. Temperature measurements at sufficiently high resolution also provide a way to estimate turbulence parameters relatively simply, building on the ideas first put forward by Osborn and Cox (1972) and Thorpe (1977).

In this context, high-resolution thermistors have been developed at the Royal Netherlands Institute for Sea Research [Nederlands Instituut voor Zeeonderzoek (NIOZ)] during the last decade. A description of the thermistors is provided in van Haren et al. (2009), even if the design has been constantly revised and improved afterward. The thermistors provide a precision that can be better than 0.5 mK (depending on the calibration). They can sample independently at frequencies up to 2 Hz for long periods of time, up to several months. These instruments have been successfully used in several oceanographic studies (e.g., van Haren and Gostiaux 2009, 2010; van Haren et al. 2012, 2015; Cimatoribus and van Haren 2015, 2016).

Negative temperature coefficient thermistors are small semiconductor sensors that exhibit large nonlinear changes in resistance as a function of temperature. It is well known that the response of thermistors to temperature is not constant (Wood et al. 1978), but rather drifts in time at a rate that depends on temperature, as well as on the age and build characteristics of the sensor. A standard procedure to compensate for this drift in an oceanographic context, however, is not well established. This is possibly because thermistors are seldom deployed for long enough periods to make this drift detectable, or because the recorded signals are most often much larger than the drift (e.g., large seasonal cycle). In most NIOZ thermistors, this drift can be detected in deep ocean observations on time scales of at least a few weeks; typical drift rates can reach up to 1 mK week−1. In this paper, we describe a physically motivated processing algorithm to compensate for this drift. While the procedure is designed and tested on NIOZ instruments alone, we argue it could be easily adapted and applied to data from other devices. The main aim of the procedure is to increase the precision of the instruments, that is, to reduce the relative drift of the temperature measured by different thermistors. While accuracy is obviously desirable as well, precision is generally the main requirement for turbulence parameter estimation, since only temperature variations (in time and in the vertical direction) are used. It will be discussed that vertical resolution is the main factor determining the applicability of the method.

In section 2, the drift of a set of NIOZ thermistors is analyzed in carefully controlled laboratory conditions. In section 3, a drift estimation algorithm for observational data is described and applied on a dataset having a very high signal-to-noise level and no missing data. The procedure is further tested using two other more challenging datasets in section 4. A summary and concluding remarks follow in section 5.

2. Laboratory analysis of the drift

The thermistors are usually calibrated at NIOZ, exploiting a custom-built calibration bath, enabling accurate calibrations in the range between −3° and +30°C. This calibration bath was originally designed to increase the temperature stepwise within this range. At each step, the temperature is maintained constant for approximately 30 min. Spatial fluctuations are below 0.1 mK in the 0.03-m-thick titanium plate at the bottom of the bath. The sensor tips, when immersed in the bath, are fixed into this titanium plate. The uniformity of the temperature in the plate has been verified during several experiments by (re)distributing the sensors randomly and by using two Sea-Bird Electronics (model 35) standard platinum high-precision thermometers. These reference thermometers are aged and calibrated by the manufacturer and have an estimated drift rate of less than 1 mK yr−1 (according to the manufacturer). The working fluid of the bath is a mix of water and glycol to prevent ice formation.

This calibration bath was adapted to analyze the drift of the sensors. A total of 40 NIOZ thermistors of different ages and designs have been used in the experiment. In particular, 10 thermistors are NIOZ3 sensors built in 2005 and 30 thermistors are more recent NIOZ4 sensors, built in 2009. The thermistors were previously stored in air at room temperature.

The drift experiment itself consisted of keeping the sensors in the calibration bath for approximately 25 days while temperature is maintained at 10.5° ± 0.1°C. Since the bath was not originally designed for these working conditions, frequent manual tuning of the bath had to be performed in order to avoid excessive departures from a constant temperature. Spatial inhomogeneity in the bath during this kind of experiment, as estimated from the two Sea-Bird thermometers, is below 2 mK. Note that this relatively high value is at least one order of magnitude larger than estimated during several calibrations performed with this apparatus. Relatively large fluctuations of temperature in time and space are due to the fact that we are using the calibration bath for a purpose for which it was not built; the bath is designed to hold temperature uniform and constant only for approximately 30 min. The results of the experiment are summarized in Fig. 1. The figure shows the temperature as measured by the thermistors and by the reference Sea-Bird thermometer (the only one available in this case and during calibration), averaged in a moving window of 2 h to remove the effect of any fast temporal fluctuation (smaller than 1 mK). The data from the first 5 days of the experiment are excluded from further analysis, due to the relatively large temperature fluctuations.

Fig. 1.

Overview of the laboratory experiment, showing the temperature measured by the Sea-Bird thermometer (black) and the temperatures measured by each thermistor (colored lines). The legend shows the unique identification number of each thermistor in the experiment, ranging from 10 to 20 for NIOZ3 sensors and from 311 to 342 for newer NIOZ4 sensors. Moving averages in 2-h-long windows are shown.

Fig. 1.

Overview of the laboratory experiment, showing the temperature measured by the Sea-Bird thermometer (black) and the temperatures measured by each thermistor (colored lines). The legend shows the unique identification number of each thermistor in the experiment, ranging from 10 to 20 for NIOZ3 sensors and from 311 to 342 for newer NIOZ4 sensors. Moving averages in 2-h-long windows are shown.

The thermistors were calibrated in the same bath after the drift estimation experiment, evaluating the response of the thermistors to changes in temperature between 5° and 14°C in steps of 1°C. This limited temperature range of the calibration bath is used to reduce the calibration time and to reduce the nonlinearity of the calibration curve, thus increasing its accuracy. Data from two NIOZ3 sensors were excluded from the analysis, having a maximum calibration error higher than 0.2 mK. This strict requirement is used to guarantee that the sensor drift can be distinguished from the calibration error during the relatively short experiment. The nonlinear response of the thermistors to temperature is fitted, using the least squares method, by the function

 
formula

where is the digitized reading of the thermistor, is the temperature measured by the Sea-Bird thermometer, and are the coefficients determined in the fit. The noise level of the instruments is below 1 mK, in most sensors by more than one order of magnitude (not shown).

Figure 2a shows the difference between the measurement of each thermistor and the reference Sea-Bird thermometer. Deviations from the reference are asymmetrically spread around zero, with more negative deviations in the first part of the experiment and more positive afterward. These systematic biases are likely due to the combined effect of the thermistor drift and of the temporal and spatial fluctuations of temperature in the bath during the long experiment. Since the calibration was performed after this experiment (after collecting the data in Fig. 2), it looks as if the drift is taking place backward in time, with a larger spread around the reference at the beginning of the experiment. The figure also shows that the signal from the thermistors has a small (≤1 mK) fluctuating component that is similar among different sensors. This shared fluctuating component is the result of small variations in time of the water temperature due to, for example, changes in the settings of the bath, opening of the bath lid to perform maintenance (spike at day 16 in Figs. 1 and 2a), and variations in room temperature, for which the bath is not able to immediately compensate.

Fig. 2.

(a) The deviations of the temperature measured by the thermistors (colors as in Fig. 1) from the temperature measured by the Sea-Bird thermometer. The black line is the shared fluctuating component discussed in the text. (b) The estimate of the thermistors drift.

Fig. 2.

(a) The deviations of the temperature measured by the thermistors (colors as in Fig. 1) from the temperature measured by the Sea-Bird thermometer. The black line is the shared fluctuating component discussed in the text. (b) The estimate of the thermistors drift.

This shared fluctuating component is estimated by taking the average of the deviation of each thermistor from the Sea-Bird thermometer. If is the temperature measured by the ith thermistor and is the reference temperature measured by the Sea-Bird thermometer, then is defined as

 
formula

where the angle brackets indicate averaging among the thermistors and the time dependence is made explicit. This shared deviation is shown as a black line in Fig. 2a.

The drift of each sensor can then be estimated as the deviation from both the Sea-Bird thermometer and this shared fluctuating component :

 
formula

The estimates of the thermistor drift based on this relation are shown in Fig. 2b. Equation (3) states that the drift can be estimated without any information on the reference temperature, but this result should obviously be taken with a grain of salt. This derives from the definition of the shared fluctuating component in Eq. (2), which implies that there is no systematic error in the temperature averaged among all the thermistors (and, equivalently, that the drift is not correlated among different sensors). If such a systematic bias were present, then Eq. (2) should include a filtering operation (e.g., time average, low pass, or other on the thermistors’ measurements), isolating such a component. Our attempts to extract such a bias, however, invariably lead to inconsistencies (>1 mK) with respect to the calibration. Since the calibrations are performed a few days after the end of this experiment, using exactly the same setup, it is highly unlikely for such inconsistencies to have a physical cause. We are thus confident that this approach is correct for this particular experimental setup.

In section 3, it will be shown that, if ocean observations are considered, only a few thermistors can be averaged to estimate the shared fluctuations and the drift. It is thus useful to provide a lower bound for the error due to the use of a limited number of thermistors. We compute the drift estimate [Eq. (3)] using, for each sensor, the average temperature of three thermistors alone, instead of all the thermistors:

 
formula

where the subscript 3 denotes that only three thermistors are used. All possible combinations (index k) of the three thermistors including thermistor i are considered.

The error is then estimated by the difference between these drift estimates and the estimate computed using the average from all the thermistors:

 
formula

The standard deviation of all the is below 0.6 mK during the whole experiment (increasing with the distance in time from the calibration). Figure 3 shows the maxima among the combinations for each thermistor, , during the whole experiment. The figure demonstrates that the drift of the sensors introduces an uncertainty in the reference estimation. It also suggests that, in order to minimize this error, the reference should be estimated at intervals of a few days. The figure also suggests that the temperature inhomogeneity estimated using the two reference thermometers may be overestimated ( is below 1 mK for all sensors at the end of the experiment).

Fig. 3.

Estimates of the error of the drift if only three thermistors are used to compute it. The colors are the same as in Fig. 1. See text for details.

Fig. 3.

Estimates of the error of the drift if only three thermistors are used to compute it. The colors are the same as in Fig. 1. See text for details.

To remove the drift estimate [Eq. (3)] from the measured temperature signal, a functional form for the drift is derived. Both Fig. 2b and previous studies (Wood et al. 1978) suggest that the thermistors, when kept at fixed temperature, approach a constant drift rate after some time. If, on the other hand, the temperature of the environment in which the thermistors are kept is changed significantly, then their drift rate changes, approaching after some time the new asymptotic value associated with the new environmental temperature. At first, the drift rate changes faster, while on a longer time scale the “equilibrium” drift rate is approached.

A few laboratory tests indicate that changes of 5°–10°C are needed to cause a detectable change in the drift rate of NIOZ thermistors. This is what happened during the experiment discussed above, when the thermistors previously kept at approximately 20°C were deployed in the bath at 10.5°C. The new asymptotic value is usually reached in a few weeks, as Fig. 2 also suggests. The change in drift rate is caused by changes in temperature rather than pressure, which is virtually the same inside the bath as outside (the response of earlier versions of the instruments had a weak dependence on pressure; see van Haren et al. 2009).

The drift rate is modeled as follows:

 
formula

with a, τ, β, and m as four constants with appropriate dimension (a, τ, and β are positive). The first term on the right-hand side of Eq. (4), m, is the constant drift rate that is approached at a constant temperature after a sufficiently long time. While the aging of sensors may actually affect m, this could not be detected in either the laboratory or field observations, likely for being too slow a process. It is probably relevant only on time scales longer than one year and for more recently built sensors. The second term on the right-hand side of Eq. (4) is a stretched exponential term. We are unaware of the previous use of this functional form in this context. Originally introduced to represent the relaxation of a capacitor by Kohlrausch (1854), the stretched exponential has possibly been most successfully used to represent mechanical relaxation in glassy materials (Phillips 1996). It represents here the relaxation toward the constant drift rate at some specific temperature.

The drift at time t is obtained by integrating Eq. (4):

 
formula

where γ is the lower incomplete gamma function , A is an amplitude factor [different from a in Eq. (4)], and is a constant bias that represents possible systematic errors. This is the functional form expected for the thermistors drift. Without loss of generality, it is fitted on the estimated drift from Eq. (3) in the form

 
formula

which, in comparison to Eq. (5), avoids the inconvenience of the explicit dependence between the two parameters of the incomplete gamma function. Of course, the two parameters can be (and in fact are) nonetheless correlated.

The fitting function is highly nonlinear, and obtaining a satisfactory fit of the data is not trivial. The fitting algorithm, defined after several tests and careful tuning, is detailed below.

  1. Let the drift of one thermistor be estimated as through Eq. (3), at a set of equally spaced discrete times , for .

  2. A least squares linear fit of the drift is performed, obtaining a first guess of the functional form of the drift , where the subscript l refers to “linear.” The coefficient of determination of this fit is .

  3. Provided that N is sufficiently large (), two fits of Eq. (6) are attempted using a simplex algorithm (Nelder and Mead 1965). The fits are performed from an initial guess with , , and . One series uses as an initial guess of , while the other series uses ; the initial guess for A is .

  4. The coefficient of determination of each of the two fits, if successful, is computed. The best fit is selected as the provisional fit of with the coefficient of determination .

  5. The amplitude A is increased by a factor of 1.5, and points 3 and 4 are iterated, always storing the best fit. The iteration is repeated for up to 40 times. The iteration stops before 40 iterations if several increases of A do not improve further.

  6. Once the iteration is completed, the nonlinear fit of is preferred over the linear one if .

This relatively elaborate fitting procedure is used to minimize the dependence on the initial conditions, by scanning a broad range of amplitudes A and two values of β. More values of β could, in principle, be used, but this increases the computational effort without any substantial improvement of the results. The final condition at step 6 is used to make sure that the nonlinear fit is a significant improvement over the linear one, explaining at least 30% of the variance not explained by the linear one.

The fitting procedure gives good results, with coefficients of determination very close to 1 with the exception of the thermistors that hardly show any drift. For those sensors, the fitting algorithm still provides reasonable results, by fitting a straight line approximately horizontal and close to the zero axis. Four examples of the results of the fit are shown in Fig. 4. The function defined in Eq. (6) provides an excellent model for the deviations from Eq. (3), even when the drift rate is rapidly changing during the experiment.

Fig. 4.

Four examples of the results of the fit described in the text, using the drift estimates shown in Fig. 2b (blue dots). The thermistor identification number is reported on top of each panel. The final results of the fitting procedure are shown as a continuous red line. If the nonlinear fit is selected as the final result, then the linear fit is shown as a dashed red line. The coefficient of determination of the linear () and, if finally selected, of the nonlinear () fit are reported in each panel. The vertical scale is different in every panel.

Fig. 4.

Four examples of the results of the fit described in the text, using the drift estimates shown in Fig. 2b (blue dots). The thermistor identification number is reported on top of each panel. The final results of the fitting procedure are shown as a continuous red line. If the nonlinear fit is selected as the final result, then the linear fit is shown as a dashed red line. The coefficient of determination of the linear () and, if finally selected, of the nonlinear () fit are reported in each panel. The vertical scale is different in every panel.

The temperature measurements from the thermistors, once compensated for this drift, deviate from the reference temperature by less than 0.1 mK. This value is likely smaller than the actual error of the overall calibration procedure. However, the fact that the thermistors measure very small fluctuations (less than 1 mK) at periods shorter than approximately one hour (not shown) suggests that inhomogeneities are indeed small in the calibration bath and that the procedure is accurate to better than 1 mK at least. A more detailed analysis of the accuracy of the calibration procedure goes beyond the scope of this article.

3. Application to ocean observational data

The laboratory analysis of the thermistor drift provides the basis for a drift compensation algorithm to be applied on ocean data. The main difficulty, when dealing with observational data, is the lack of a reference temperature measurement. In other words, the term in Eq. (3) has to be substituted. This average among the thermistors cannot be used since thermistors, usually deployed at different depths, measure on average different temperatures. We illustrate the approach using the dataset described briefly in Table 1 (mooring 1) and in more detail in Cimatoribus and van Haren (2015). During this deployment, the thermistors worked particularly well, with very low noise levels, accurate calibrations, and no missing data.

Table 1.

Description of the moorings from which the data used in this work were collected. Acronym h.a.b. is the height above the bottom.

Description of the moorings from which the data used in this work were collected. Acronym h.a.b. is the height above the bottom.
Description of the moorings from which the data used in this work were collected. Acronym h.a.b. is the height above the bottom.

a. First guess of the reference temperature profile

The drift identification procedure is composed of two steps. During the first step, measurements from each thermistor are divided and averaged over consecutive, independent time windows. The length of the time window can vary depending on the oceanographic conditions. The guiding principle is that the window should be as short as possible, averaging only the isothermal displacements due to gravity waves and unstable turbulent overturns.1 The drift rate changes quickly after the deployment, and shorter windows better resolve these variations. In practice, based on tests using different datasets and configurations, and after comparison with nearby shipborne conductivity–temperature–depth (CTD) profiles, we find that a length of one day generally gives the best results. Longer windows provide comparable results if the first weeks after the deployment are not considered.

Figure 5 shows the measurements from the thermistors (mooring 1 of Table 1), averaged over one day (9 July 2013, randomly chosen). The main figure shows the entire depth range of the mooring, while the inset is a zoom-in of a smaller interval. The black line shows the measurements from the thermistors. Differently from the laboratory experiment described in the previous section, some thermistors use a calibration performed before deployment, some use a calibration performed after deployment, and some use a combination of the two. Multiple calibrations are combined, averaging the coefficients of Eq. (1). This is a common situation, since the bath cannot hold all the sensors at the same time. Time constraints are often a limiting factor as well.

Fig. 5.

The black line shows the temperature measured by the thermistors of mooring 1 (see Table 1), averaged over one day (9 Jul 2013). The gray squares mark thermistors that are excluded from the fit of the polynomial. Blue crosses show the smoothed temperature profile, and the red line is the final estimate of temperature after the full drift compensation procedure is applied. See text for further details. The inset is an enlarged view over a short depth interval. The acronym h.a.b. stands for height above the bottom.

Fig. 5.

The black line shows the temperature measured by the thermistors of mooring 1 (see Table 1), averaged over one day (9 Jul 2013). The gray squares mark thermistors that are excluded from the fit of the polynomial. Blue crosses show the smoothed temperature profile, and the red line is the final estimate of temperature after the full drift compensation procedure is applied. See text for further details. The inset is an enlarged view over a short depth interval. The acronym h.a.b. stands for height above the bottom.

The temperature profile in Fig. 5 (black line) clearly shows that some thermistors have suffered from a strong drift, and overall the profile is not as smooth as we would expect from a 1-day average computed from 1-Hz data. To make progress, we have to assume that, on average, the overall stratification is still captured by the thermistors, and we must make a guess of this average stratification within each averaging window. Various interpolation methods have been tested, including polynomial fits of different order, piecewise cubic Hermite interpolating polynomials, Fourier filtering, and smoothed splines. Overall, we find that a polynomial fit, generally of degree 3 or 4, gives a good estimate of the average profile with a minimum amount of configuration (the degree of the polynomial alone). Also smoothed splines give good results, better in some ocean stratification conditions, but they require more careful tuning. The smoothed profile of this example, from a least squares fit of a polynomial of degree 4, is shown in the figure with blue crosses.

Extra care has to be taken when, for some reason, the data from one or more thermistors are missing, or when some thermistors have a systematic bias with respect to the others (in Fig. 5, one such biased thermistor is visible at approximately 35 m above the seafloor). In this latter case, we found it useful to exclude these thermistors from the fit (thermistors marked by a gray square in Fig. 5), to prevent biases of the estimated average stratification. A practical way to do this is setting a threshold on the value of the second derivative of temperature with respect to depth, estimated from the nonsmoothed profile. In Fig. 5, a threshold of 0.08°C m−2 was used.

A first guess of the drift of a thermistor within a single averaging time window is the difference between the thermistor measurement and the smoothed profile. This first guess of the drift, of thermistor i in the time window j, is represented by .

b. Identification of the time-dependent drift

Compensating the drift by subtracting from a thermistor measurement may be acceptable if accuracy is not a concern and if a rather short portion of the data is used (comparable to the averaging window length). In principle, we could try to extract the thermistor drift by fitting the function [Eq. (6)] directly on the sequence of of each thermistor. However, this is generally unfeasible, since the drift signal is small compared to the fluctuations due to natural variability (see section 4 for an example of the opposite).

A more reliable estimation (and compensation) of the drift is possible by following an approach inspired by the one used in the laboratory. Figure 6a shows of three nearby thermistors in the mooring, respectively 21.1 (red), 21.8 (blue), and 22.5 m (green) above the seafloor. These three thermistors can also be identified in the inset of Fig. 5, with the thermistor 21.8 m above the seafloor positively biased with respect to the others. As in the laboratory, we recognize that there is a similar fluctuating component among the three nearby thermistors. In this case, these fluctuations are the effect of the variability of stratification due to currents and waves. A first guess of this shared fluctuating component is computed as

 
formula

where the overline indicates temporal averaging. In Eq. (7), data from three thermistors are used, , i, and —that is, the three nearby thermistors centered on the one we are studying. The subtraction of the time-averaged removes constant systematic biases (e.g., due to calibration) from the shared fluctuating component. The second guess of the drift is then given by . This second guess of the drift is then fitted using the same procedure used for the laboratory data and described in section 2. We indicate the function resulting from the fit as [either a straight line or the function of Eq. (6)]. We use this preliminary fit for detrending the sequence:

 
formula

the detrended are averaged in groups of three, similarly to what is done in Eq. (7), to obtain the shared fluctuating component :

 
formula

The detrending is needed to compensate for the very different magnitude that the drift rate can take in different sensors and systematic biases that cannot be averaged out among only three sensors. This issue was not relevant in the laboratory, since averaging could be performed among several sensors therein, minimizing the influence of individual, above average, drift rates.

Fig. 6.

Summary of the procedure for the identification and fit of the thermistor drift in observational data. (a) The (above the seafloor) for three nearby thermistors: 21.1 m (red), 21.8 m (blue), and 22.5 m (green) lines. The is shown as a thick black line. (b) The for the thermistor 21.8 m above the seafloor are shown again as blue crosses [corresponding to the blue line in (a)], and is shown as black circles. Note that the vertical scale is different in the two panels. Two fits of the drift are shown in red, a linear one (dashed line) and one using the function defined in Eq. (6) (continuous line). The coefficients of determination for the two fits are reported.

Fig. 6.

Summary of the procedure for the identification and fit of the thermistor drift in observational data. (a) The (above the seafloor) for three nearby thermistors: 21.1 m (red), 21.8 m (blue), and 22.5 m (green) lines. The is shown as a thick black line. (b) The for the thermistor 21.8 m above the seafloor are shown again as blue crosses [corresponding to the blue line in (a)], and is shown as black circles. Note that the vertical scale is different in the two panels. Two fits of the drift are shown in red, a linear one (dashed line) and one using the function defined in Eq. (6) (continuous line). The coefficients of determination for the two fits are reported.

This relatively complex procedure guarantees a good estimate of the drift also at the beginning of a dataset, when the drift rate is relaxing to its new value set by the temperature of the deployment environment. The shared fluctuating component of the thermistor at 21.8 m is shown as a black line in Fig. 6a, where we can see that the detrending successfully avoided the inclusion of the drift in the fluctuating component. By removing this shared fluctuating component from the initial sequence of , reported in Fig. 6b with blue crosses, we obtain the final estimate of the drift, shown as black circles in the same panel:

 
formula

This estimate of the drift can finally be fitted using again the same procedure, obtaining , the fitted functional form of the drift. The linear (dashed red line) and stretched exponential (continuous red line, the one actually used) fits are shown in Fig. 6b, with their respective coefficients of determination.

The thermistors drift compensation is performed by repeating this procedure for each thermistor and subtracting the results of the fit from the temperature measurements. For the top and bottom thermistors, only the nearest sensor is used to identify the shared fluctuating part. The good results shown in Fig. 6 are in no way exceptional but are instead representative of what is obtained with the other thermistors in the mooring. In Fig. 5, the red line shows the temperature profile after the complete drift compensation procedure (1-day average). Note that the drift compensation procedure retains finer structure in the stratification, in comparison to the first guess of the profile.

The overall effect of the drift compensation can be qualitatively evaluated from Fig. 7. Figures 7a and 7b show temperature measurements, and Figs. 7c and 7d show the vertical temperature gradient during one hour on the first day of deployment of the instruments. The full-resolution data are plotted without any interpolation or smoothing. This time interval is chosen because, while the value of the drift is larger at the end of the deployment, the drift compensation is more challenging in the first part of the dataset, when the drift rate is usually larger and changing in time (adapting to the ocean conditions after deployment). In Figures 7a and 7c, the differential drift of some thermistors is revealed by several horizontal bands, more visible in the parts of the mooring having weaker stratification and in the vertical gradient. The drift compensation procedure is very effective in removing virtually any trace of these bands, without affecting the overall stratification or patterns. Further comparisons of the measurements with independent CTD profiles, collected in the same area after recovering the mooring, confirm that the average stratification measured by the thermistors is consistent with that measured by the CTD. The temperature stratification estimated by the thermistors, after drift compensation, is 2.4 mK m−1 (average during the day of the CTD survey), while the one from the CTD is 2.7 mK m−1 (average of one CTD profile in the same depth interval as the moored thermistors, nearby location). The gradient is highly variable throughout the water column, with a standard deviation of 5 mK m−1 in both estimates.

Fig. 7.

Temperature measurements during one hour, at mooring 1 of Table 1, from 0000 to 0100 UTC 14 Apr 2013. The temperature measurement (a) before and (b) after the drift compensation procedure. The vertical temperature gradient (c) before and (d) after the drift compensation procedure. On the vertical axes, three arrows mark the depths of the three thermistors used in Fig. 6a, with corresponding colors. The acronym h.a.b. stands for height above the bottom.

Fig. 7.

Temperature measurements during one hour, at mooring 1 of Table 1, from 0000 to 0100 UTC 14 Apr 2013. The temperature measurement (a) before and (b) after the drift compensation procedure. The vertical temperature gradient (c) before and (d) after the drift compensation procedure. On the vertical axes, three arrows mark the depths of the three thermistors used in Fig. 6a, with corresponding colors. The acronym h.a.b. stands for height above the bottom.

Overall, the essential constraint of the method is that the smoothing does not introduce a systematic bias. How to prevent such biases depends on the particular dataset analyzed, but as long as the thermistor spacing is small compared to the isopycnal displacement, it is reasonable to assume that any fine structure will be smoothed out by averaging in a sufficiently long time window. While the smoothing procedure may not always correctly identify the mean stratification within every window, the procedure outlined above is able to identify these biases by considering multiple thermistors together, and the fit of Eq. (6) further imposes a physically based constraint on the estimate of the drift.

We attempt a more quantitative analysis of the effect of changing the vertical resolution by considering again the drift of one sensor—the same one considered in Fig. 6—21.8 m above the seafloor. We compute the drift estimate for this thermistor using increasingly distant sensors, rather than the two nearest ones alone. To be more specific, the sum in Eqs. (7) and (8) is substituted with with . Every estimate at larger m mimics a lower vertical-resolution measurement. The results of this exercise are shown in Fig. 8, which reports the deviation from the original, full-resolution (, 0.7 m), drift estimate for (with a vertical resolution of 1.4, 2.1, 2.8, and 4.2 m, respectively).

Fig. 8.

Estimate of the error on the drift reducing vertical resolution. Different colors refer to different vertical resolutions (see legend), obtained by subsampling in the vertical direction the dataset of mooring 1 (Table 1). The drift estimate is computed for the sensor 21.8 m above the seafloor; see text for details.

Fig. 8.

Estimate of the error on the drift reducing vertical resolution. Different colors refer to different vertical resolutions (see legend), obtained by subsampling in the vertical direction the dataset of mooring 1 (Table 1). The drift estimate is computed for the sensor 21.8 m above the seafloor; see text for details.

From the figure we conclude that the method produces consistent results (with a precision of approximately 1 mK) for vertical spacing up to approximately 3 m. For larger spacing, deviations from the full-resolution estimate are larger, and large fluctuations and a decreasing trend in time are evident. Only the 4.2-m resolution is shown as an example, but larger spacings provide similar or worse results.

While this test gives some guidance for estimating the vertical-resolution requirements in similar conditions, we stress that these results should not be taken as general. The resolution requirements will change markedly with—to name a few examples—stratification, coherence of vertical motions throughout the mooring, presence of strong horizontal gradients, etc. Ultimately, only a careful analysis of the results obtained from a particular mooring can provide an answer.

4. Examples of applications

We now briefly present two other applications of the described procedure, using data from two other moorings. First, we consider data from mooring 2 in Table 1, deployed immediately after mooring 1 at a nearby, deeper location. This dataset was used and described in more detail in Cimatoribus et al. (2014) and van Haren et al. (2015). This mooring used 140 thermistors of different ages, attached with a 1-m spacing between them. Because of the deeper location, stratification is weaker. Furthermore, battery failures and calibration issues afflicted this deployment, with 7 out of 140 thermistors providing either unreadable or unusable data (white stripes in Figs. 10a and 10b). These instrumental issues, in combination with the weaker stratification, make this a more challenging test for the drift compensation procedure. The drift compensation procedure is performed in the same way as for mooring 1, and it performs time averaging in windows of 1 day but uses a polynomial of degree 3 (instead of 4) as a guess of the real stratification, since stratification is closer to linear at this location. The drift of thermistors next to a missing one is estimated using the available data alone; that is, no interpolation is performed, and no sensors other than the nearest ones are used.

Figure 9 shows an example of the drift estimation and fit for the thermistor 51 m above the seafloor. We see that also in this case the drift is extracted and fitted reliably, with the negative drift rate identified from the initial estimates (blue line in Fig. 9a) despite a nearby thermistor having a stronger, opposite trend (red line in Fig. 9a). Figure 10 confirms that the drift compensation is effectively removing the drift from the observations, and it shows that the procedure is not significantly affected by the lack of data from some of the thermistors.

Fig. 9.

As in Fig. 6, but for the dataset from mooring 2 in Table 1. (a) Data from thermistors at a height above the seafloor of 55 (red), 56 (blue) and 57 m (green). (b) Drift estimation for thermistor 56.

Fig. 9.

As in Fig. 6, but for the dataset from mooring 2 in Table 1. (a) Data from thermistors at a height above the seafloor of 55 (red), 56 (blue) and 57 m (green). (b) Drift estimation for thermistor 56.

Fig. 10.

As in Fig. 7, but for the dataset from mooring 2 in Table 1. Data from 1200 to 1300 UTC 13 Aug 2013. The acronym h.a.b. stands for height above the bottom.

Fig. 10.

As in Fig. 7, but for the dataset from mooring 2 in Table 1. Data from 1200 to 1300 UTC 13 Aug 2013. The acronym h.a.b. stands for height above the bottom.

The last dataset considered has been collected between years 2006 and 2007 using older, NIOZ3, thermistors. This dataset was described in more detail in van Haren et al. (2009), van Haren and Gostiaux (2009), and Gostiaux and van Haren (2012). The mooring was anchored in deep water (5275 m), but the thermistors were attached below a top buoy at a shallower depth, approximately 1500 m. The water column is more stratified at this depth than in the other two moorings. However, the calibration had to be performed in the field at the time of the deployment, using concurrent CTD measurements, since the calibration bath had not been built yet. This posed stronger limitations on the accuracy of the thermistor measurements. Furthermore, some of these older thermistors suffered from an extreme drift (in particular, the thermistor at 1410 m in Figs. 11 and 13a). Finally, during the long deployment, some thermistors failed due to battery depletion, so that not every time series has the same length.

Fig. 11.

As in Fig. 5, but for the dataset from mooring 3 in Table 1. Data are 1-day averages from 27 Oct 2006.

Fig. 11.

As in Fig. 5, but for the dataset from mooring 3 in Table 1. Data are 1-day averages from 27 Oct 2006.

Despite these issues, the drift-compensating procedure can still be applied successfully, with minimal changes. In particular, a smoothed spline is used in this case rather than a polynomial fit for the initial estimation of stratification, since the vertical temperature gradient at this depth has more vertical structure than in the deeper moorings (see Fig. 11). The drift estimation is then performed first on the thermistor with the strong drift at 1410 m alone. This thermistor has a drift of more than 1°C in one year, and for this reason its drift can be estimated without removing the fluctuations, which are much smaller than the total drift. Then, the drift of all other thermistors can be estimated as done for the other moorings. This two-step procedure enables the correct estimation of the drift of the sensors adjacent to the one with extreme drift, which would otherwise suffer from a biased estimation of the fluctuating part.

Figure 12 shows the drift identification and fit for the thermistor at 1450 m. Nearby thermistors (red and green lines) stop measuring before this thermistor (black line), but this does not affect the drift estimation significantly, because the drift rate has approximately reached its asymptotic value by the time the other thermistors stop. Note that the fluctuations (black line) show a clear positive trend, which can be attributed to the passage of a warm mesoscale structure at the mooring location when the complete dataset is considered (van Haren et al. 2009). This trend is correctly included in the common fluctuations and thus removed from the drift estimate (i.e., it is left in the natural temperature variation signal).

Fig. 12.

As in Fig. 6, but for the dataset from mooring 3 in Table 1. (a) Data from thermistors at a depth of 1452 (red), 1449.5 (blue), and 1447 m (green).

Fig. 12.

As in Fig. 6, but for the dataset from mooring 3 in Table 1. (a) Data from thermistors at a depth of 1452 (red), 1449.5 (blue), and 1447 m (green).

The fit represents the drift very well also in this case. We note that these older thermistors generally have higher drift rates and that the drift is often fitted sufficiently well by a straight line in this dataset. Despite the challenges encountered in processing this last dataset, the results of the drift compensation are convincing (Fig. 13), even at the beginning of the deployment. The temperature stratification is negative in large parts of this dataset due to the density compensation of salinity variations at this location.

Fig. 13.

As in Fig. 7, but for the dataset from mooring 3 in Table 1. Data from 0000 to 0100 UTC 11 Jun 2006.

Fig. 13.

As in Fig. 7, but for the dataset from mooring 3 in Table 1. Data from 0000 to 0100 UTC 11 Jun 2006.

5. Summary and concluding remarks

We presented a characterization of the drift of the temperature sensors designed at NIOZ and used in several oceanographic moorings. The characterization is done in the custom-built calibration bath of NIOZ. The laboratory experiment supports the hypothesis that the drift rate relaxes in time toward an asymptotic value, which is itself a function of temperature. Using this information, we develop a procedure to estimate the drift rate of thermistors deployed in the field, where no reference temperature measurements are available. This procedure implies a guess of the average stratification computed from the thermistor measurements in independent time windows, by means of a smoothed interpolation. The drift estimated from this first guess is then fitted on the basis of the expected relaxation toward an asymptotic drift rate. To successfully perform this fit, fluctuations due to natural variability of stratification have to be removed, after estimating them by exploiting the combined information from nearby sensors.

This procedure has been described emphasizing the application on data collected using NIOZ thermistors. However, it is readily applicable to other cases where temperature is recorded by multiple thermistors at different depths. This is true in particular because the functional form of the drift seems to be similar in different instruments. The limiting factor is vertical resolution, which must be sufficiently high to guarantee that fluctuations due to waves and other natural variability can be identified among different sensors.

Acknowledgments

The authors thank the technicians at NIOZ, especially Martin Laan, for the technical support during all measurements, particularly concerning the configuration, running, and tuning of the calibration bath. The authors would also like to thank the crew of the R/V Pelagia for enabling the collection of the oceanographic measurements.

REFERENCES

REFERENCES
Cimatoribus
,
A. A.
, and
H.
van Haren
,
2015
:
Temperature statistics above a deep-ocean sloping boundary
.
J. Fluid Mech.
,
775
,
415
435
, doi:.
Cimatoribus
,
A. A.
, and
H.
van Haren
,
2016
:
Estimates of the temperature flux–temperature gradient relation above a sea floor
.
J. Fluid Mech.
,
793
,
504
523
, doi:.
Cimatoribus
,
A. A.
,
H.
van Haren
, and
L.
Gostiaux
,
2014
:
Comparison of Ellison and Thorpe scales from Eulerian ocean temperature observations
.
J. Geophys. Res. Oceans
,
119
,
7047
7065
, doi:.
Gostiaux
,
L.
, and
H.
van Haren
,
2012
:
Fine-structure contamination by internal waves in the Canary Basin
.
J. Geophys. Res.
,
117
,
C11003
, doi:.
Kohlrausch
,
R.
,
1854
:
Theorie des elektrischen Rückstandes in der Leidener Flasche
.
Ann. Phys.
,
167
,
179
214
, doi:.
Nelder
,
J. A.
, and
R.
Mead
,
1965
:
A simplex method for function minimization
.
Comput. J.
,
7
,
308
313
, doi:.
Osborn
,
T. R.
, and
C. S.
Cox
,
1972
:
Oceanic fine structure
.
Geophys. Fluid Dyn.
,
3
,
321
345
, doi:.
Phillips
,
J. C.
,
1996
:
Stretched exponential relaxation in molecular and electronic glasses
.
Rep. Prog. Phys.
,
59
,
1133
, doi:.
Thorpe
,
S. A.
,
1977
:
Turbulence and mixing in a Scottish loch
.
Philos. Trans. Roy. Soc. London
,
A286
,
125
181
, doi:.
van Haren
,
H.
, and
L.
Gostiaux
,
2009
:
High-resolution open-ocean temperature spectra
.
J. Geophys. Res.
,
114
,
C05005
, doi:.
van Haren
,
H.
, and
L.
Gostiaux
,
2010
:
A deep-ocean Kelvin-Helmholtz billow train
.
Geophys. Res. Lett.
,
37
,
L03605
, doi:.
van Haren
,
H.
,
M.
Laan
,
D.-J.
Buijsman
,
L.
Gostiaux
,
M. G.
Smit
, and
E.
Keijzer
,
2009
:
NIOZ3: Independent temperature sensors sampling yearlong data at a rate of 1 Hz
.
IEEE J. Oceanic Eng.
,
34
,
315
322
, doi:.
van Haren
,
H.
,
L.
Gostiaux
,
M.
Laan
,
M.
van Haren
,
E.
van Haren
, and
L. J. A.
Gerringa
,
2012
:
Internal wave turbulence near a texel beach
.
PLoS One
,
7
,
e32535
, doi:.
van Haren
,
H.
,
A.
Cimatoribus
, and
L.
Gostiaux
,
2015
:
Where large deep-ocean waves break
.
Geophys. Res. Lett.
,
42
,
2351
2357
, doi:.
Wood
,
S. D.
,
B. W.
Mangum
,
J. J.
Filliben
, and
S. B.
Tillet
,
1978
:
An investigation of stability of thermistors
.
J. Res. Natl. Bur. Stand.
,
83
,
247
263
, doi:.

Footnotes

Publisher’s Note: This article was revised on 26 July 2016 to correct the incomplete inclusion of the first author’s current affiliation to the title page of the article.

a

Current affiliation: Ecological Engineering Laboratory, École Polytechnique Fédérale de Lausanne, Lausanne, Switzerland.

1

If this first guess of the drift is, for some reason, the only one that can be computed, then it is more reasonable to use a longer time window, in order to average any internal wave motion at inertial and tidal frequencies.