Methods for measuring waves and winds from a Wave Glider autonomous surface vehicle (ASV) are described and evaluated. The wave method utilizes the frequency spectra of orbital velocities measured by GPS, and the wind stress method utilizes the frequency spectra of turbulent wind fluctuations measured by an ultrasonic anemometer. Both methods evaluate contaminations from vehicle motion. The methods were evaluated with 68 days of data over a full range of open ocean conditions, in which wave heights varied from 1 to 8 m and wind speeds varied from 1 to 17 m s−1. Reference data were collected using additional sensors on board the vehicle. For the waves method, several additional datasets are included that use independently moored Datawell Waverider buoys as reference data. Bulk wave parameters are determined within 5% error with biases of less than 5%. Wind stress is determined within 4% error with 1% bias. Wave directional spectra also compare well, although the Wave Glider results have more spread at low frequencies.
The Wave Glider is an autonomous surface vehicle that uses wave motion for propulsion. A surface float is connected by a tether, typically 8 m in length, to a subsurface body with a series of six wings. With the passage of each wave, the subsurface float glides forward but not backward because the wings feather to prevent it. A rudder on the subsurface body controls navigation. The subsurface body tows the surface float, such that the forward motion is in any desired direction (i.e., propulsion is achieved both upwave and downwave, or across-wave). The surface float contains the main electronics for telemetry, command/control, and scientific instrumentation. When prepared for deployment, as in Fig. 1, the surface float and subsurface body are bundled together; a mechanical release separates them once in the water. The Wave Glider is commercially produced by Liquid Robotics, Inc., based in Sunnyvale, California. The current version is the SV3, of which over 150 units have been produced for a range of applications, from oceanographic surveying to port security.
The following work describes and evaluates methods to make scientific-grade measurements of waves and winds from the Wave Glider. For waves, the challenge is that the surface float is not truly surface following, at least not in terms of the pitch–roll–heave signals used by traditional wave buoys, and also that the vehicle is extracting propulsion from the waves. The solution herein is to use global positioning system (GPS) measurements of horizontal motion as the raw wave data signal, instead of pitch–roll–heave, in a method termed GPSwaves (Herbers et al. 2012). This method is agnostic to the orientation of the surface float as well as the details of the vehicle response to surface slopes. For wind stress, the challenge also is motion contamination along with potential biases of measuring at low heights above the water, such that the measurements may be within the wave-coherent boundary layer (e.g., Grare et al. 2013; Edson et al. 2013; Hara and Sullivan 2015). The solution herein is the use a fast sampling sonic anemometer (Gill Instruments WindMaster) to observe the inertial subrange of wind turbulence at frequencies greater than those of motion contamination (Edson et al. 1991).
These methods are assessed in a series of evaluation datasets, using the commercially produced Datawell wave sensor (Datawell BV, the Netherlands) and Airmar ultrasonic anemometer (Airmar Technology Corp., Milford, New Hampshire) as reference data. The wind datasets are limited, relative to the wave datasets, and model reanalysis products are used to bolster the validation of this portion.
There have been several previous investigations of wave and wind measurements from the Wave Glider, particularly using the earlier SV2 model. Focusing on waves only, Wang and Allard (2012) find that the measurements are in good agreement with nearby National Data Buoy Center (NDBC) measurements. They note that wave spectra measured with the Wave Glider are slightly broader with some bias at low frequencies. More recently, Lenain and Melville (2014) used a Datawell MOSE-G sensor mounted on board a Wave Glider to measure wave heights up to 10 m and winds up to 37 m s−1 in Tropical Cyclone Freda (2012). These measurements agreed well with hindcast models of the storm conditions with adjustment of the winds to the standard 10-m reference height. Lenain and Melville (2014) also include an appendix comparing results from the on board Datawell unit with an independently moored Datawell Waverider buoy from the Coastal Data Information Program (CDIP). Although the comparison dataset is quite limited (2 days with a range of wave heights from 0.7 to 1.4 m), Lenain and Melville (2014) conclude that the Wave Glider measurements are in good agreement with the reference. Very recently, Mitarai and McWilliams (2016) presented measurements of winds up to 32 m s−1 during Typhoon Danas and suggested that the 10-min averaged wind speeds were in qualitative agreement with model reanalysis winds, despite the known complexities of the atmospheric boundary layer in the presence of surface waves (e.g., Hara and Sullivan 2015).
The present study expands upon these previous efforts using larger evaluation datasets covering full ocean conditions. The present work also describes and evaluates newly implemented methods that run on board the Wave Glider as directly integrated sensors and software packages (as opposed to third-party sensors with their own processing, such as the Datawell unit or the Airmar weather station). The wave spectral quantities and bulk wave parameters are compared using standard skill metrics (O’Reilly et al. 1996; Krogstad et al. 1999; Collins 2012). The wind stress results are compared to bulk parameterizations (i.e., Smith 1980; Fairall et al. 2003). Details of wave directional spectra will be difficult to compare (Collins 2012), even with standards such as the Wave Eval Tool (WET; Coastal Data Information Program 2009), but specific examples will demonstrate the more subtle results. It is worth noting at the outset that differences of 5%–10% in bulk wave parameters are common, even from sensors on the same platform (Bender et al. 2010).
Upon estimating waves and winds, a final evaluation and synthesis of the results will apply the wind-wave equilibrium hypothesis of Phillips (1985). This framework is one of the motivations for measuring winds and waves simultaneously on a single platform. The equilibrium hypothesis describes the short wave scales (i.e., high frequencies) that are directly forced by the wind stress. These waves scales are in a balance between 1) wind input, 2) nonlinear energy transfers within the wave spectrum, and 3) wave dissipation by whitecapping. The wave energy spectrum in this range is predicted by
where g is gravity and is the wind friction velocity, from which the wind stress is . Thus, a measure of can be used to infer independently of a direct wind measurement (Juszko et al. 1995; Thomson et al. 2013). The remaining parameters are and . The actual value of is a function of the relative angle between the waves and the wind, which will be calculated following Thomson et al. (2013). Although equilibrium is commonly observed in the high-frequency tail of wave spectra, there are notable deviations and confounding effects of swell that are topics of active research (e.g., Kukulka et al. 2013; Lenain and Melville 2017; Takagaki et al. 2012). Progress in this area will require simultaneous measurements of surface waves and wind stress, as demonstrated herein.
a. Wave measurements and processing
The approach to measuring waves from the Wave Glider follows closely the work of Herbers et al. (2012), as applied to other nonspherical surface objects, such as the Surface Wave Instrument Float with Tracking (SWIFT) drifter (Thomson 2012). The primary raw data for this method are GPS-measured horizontal velocity components: u (east) and υ (north). By using the phase of the GPS carrier wave to determine a Doppler-shifted velocity, the latest generation of GPS receivers reports velocities with a precision of a few centimeters per second (cm s−1). These velocities are used to measure the wave orbital motions. The approach is agnostic to the orientation (pitch, roll, heading) of the surface float, and it requires only that the float move laterally with the orbits of the waves. Supplemental inertial motion unit (IMU) data are include to help resolve wave direction, by calculating the phase of vertical acceleration a relative to the lateral velocity components u and υ.
1) Raw wave data collection
The raw wave data are collected using a Microstrain 3DM GX3-35 sensor mounted on a mast 0.65 m above the Wave Glider float. The GX3 sensor integrates GPS and IMU data into a single serial output. No postprocessing or real-time corrections [i.e., postprocessing kinematic (PPK) or real-time kinematic (RTK)] are made on the raw GPS data. Data are collected at a sampling frequency Hz for bursts lasting 30 min and then processed on board the vehicle using a C++ library running on the sensor management computer (SMC). The duration of the burst can be changed as an input to the onboard routines; 30 min is a recommended minimum. The raw data are written to binary files and then the select variables of are passed to the onboard processing.
2) Raw wave data preprocessing
The raw are prepared for spectral analysis by first high-pass filtering in the time domain using the classic filter
which is applied, for example, to the velocity component u as follows:
This also is applied to υ and a. The time between observations is , and RC is a time constant for the filter, such that each observation is filtered based on the previous observation . The filter is intended to avoid leakage in low frequencies caused by vehicle navigation and drift, and this resistance–capacitance (RC) filter cutoff can also be set as an input to the onboard routines. The actual cutoff frequency is , and is the default value that filters signals at frequencies lower than 0.0455 Hz (i.e., wave periods longer than 22 s).
Next, the raw data are despiked by removing any points that are more than N standard deviations from the mean, where the N is an input to the routines, with a default setting of . The RC filter and despiking do not affect the majority of measurements, but they help prevent spurious large values in cases with very low signal amplitude (i.e., small waves).
3) Spectral processing
The wave data are then parsed into windows of length with 50% overlap. These windows are detrended (linear), tapered (Hanning), and rescaled to preserve the original variance in each window. The fast Fourier transform (FFT) of each variable is calculated to produce complex Fourier coefficients at each frequency f:
where the frequency range is from the lowest frequency of to the Nyquist frequency of (the zero frequency, which corresponds to the mean, is removed). These coefficients are then combined to determine the autospectra and the cross-spectra of the signals,
where the asterisk () indicates the complex conjugate, and the frequency dependence has been dropped from the notation. The autospectra are real, while the cross-spectra are complex. The factor of 2 in the numerator arises to preserve variance when using only one side of a double-sided FFT. The factor in the denominator creates spectral densities, using the number of points in the FFT and the sampling frequency (the raw FFTs were not normalized by the number of points in the previous equations). The variance of the raw signal is distributed across a range of frequencies spanning , such that integrating the autospectra recovers the original variance of the time series (i.e., Parseval’s theorem). Finally, auto- and cross-spectra from the 128-s windows are ensemble averaged to form spectral densities for each 30-min burst that have approximately 28 degrees of freedom.
4) Spectral wave energy
Scalar wave energy density spectra are computed purely from the autospectra of the orbital velocities, using linear wave theory, where velocities relate to wave amplitude A by the relation in deep water (i.e., circular orbits). In shallow water there is a correction for the ellipticity of the orbits but that is ignored here, as most Wave Glider operations are in deep water. The wave energy spectrum, which is proportional to at each frequency, is thus
This conversion to elevation variance from velocity variance amplifies noise at low frequencies, with a known problem when measuring small amplitude low-frequency waves (with insufficient signal-to-noise ratio).
For efficiency of telemetry, only the first 49 frequency bands of this result are retained, and thus the frequencies reported are to Hz (which correspond to wave periods from 85 to 2.5 s) at a frequency resolution of 0.0078 Hz. The results are trimmed later to exclude infragravity frequencies ( Hz, or wave periods longer than 20 s).
5) Antenna height corrections
The resulting wave energy spectra are corrected for antenna height, where the height h is an input to the onboard routines (and setting would apply no correction). The correction is made by removing energy from each frequency according to the amount of lateral displacement that would be added by a tilting antenna above the surface. This extra lateral displacement would be interpreted as an additional vertical displacement, because of the assumption of circular orbits. Using the spectral mean square slope and the small angle approximation, the horizontal projection of a tilted antenna at h can be removed at each frequency of E using
This correction is generally small. For the simplified case of waves at Stokes limiting steepness of 0.443 and the Wave Glider antenna height of m, the correction is 0.08 m2 Hz−1. Given that peak wave energy densities are often in the range of 1–10 m2 Hz−1, this correction is almost negligible. Note this antenna height (with tilting) correction is much less important than it is for accelerometer-based wave measurements (e.g., Bender et al. 2010).
6) Spectral wave directions
Directional information is obtained from the cross spectra, in particular the phasing of the signals (which is related to the real vs imaginary parts of the cross spectra). Following the methods of Kuik et al. (1988) and the adaptations similar to Herbers et al. (2012), the directional moments are obtained for each frequency. These are normalized quantities, and here it is really the phasing of the and signals that matters, not the magnitude. The directional moments are used to calculate a dominant direction and a spread for each frequency band. The directional moments also can used to reconstruct a full directional spectrum with model estimators, such as the maximum entropy method (MEM; Lygre and Krogstad 1986). As the definitions and usage of are well documented in the literature, the specifics relevant to the Wave Glider are given in the appendix.
7) Quality control of spectral results
The scalar energy spectra obtained by these methods are often contaminated by spurious energy at swell frequencies, nominally at Hz. This likely is the result of the conversion from velocity spectra to heave spectra [Eq. (9)], as well as slow drift and navigation signals not removed by the RC filter or the detrending. There may also be leakage into the swell bands from the infragravity bands ( Hz), which are known to have high energy levels when sampled from quasi-Lagrangian platforms (Herbers and Janssen 2016). To avoid these effects, swell frequencies with large directional spreads are rejected in a quality-control step of the onboard processing routines. The default limit is in the criterion , though this spread cutoff can be adjusted as an input to the onboard routines. This limit can also be applied using (see the appendix). The energy densities at frequencies that fail this quality criteria are replaced by a simple linear interpolation between neighboring frequencies passing the criteria, including the assumption of zero energy density at frequencies lower than typical swell (i.e., at Hz). No more than 5 frequencies bands, out of 49, are replaced by nonzero interpolation.
8) Doppler correction for vehicle navigation
Ocean wave measurements from moving platforms typically need to be corrected for Doppler shifting of the wave signal (e.g., Collins et al. 2017). The correction is an adjustment between the intrinsic frequency, which is given by dispersion, and the absolute frequency, which is conserved as follows (in deep water):
Here is the vector wavenumber and is the average vector velocity of the platform. For opposing platform velocities, waves pass more frequently. For following platform velocities, waves pass less frequently. This redistributes energy in the frequency spectrum slightly, but it does not change the total energy (or the bulk wave height). A more complicated effect to consider would be changes to the vehicle dynamic response with different speeds , but this is beyond the scope of the present work.
Implementing a Doppler correction on board the GPSwaves routine would be an iterative process, as was done in Zippel and Thomson (2017), wherein the spectral wave directions would be estimated to get and then the corrected frequency would be applied in Eq. (9) until convergence within some tolerance was reached. This step is not included in the present version of GPSwaves, as the corrections would be only a few percent for the typical vehicle speeds. Furthermore, the average vehicle speed over a 30-min burst is often poorly defined, because the vehicle may navigate (turn, etc.) during the burst. The Wave Glider tends to advance more slowly when opposing waves, which would reduce corrections in the foreshortening cases relative to the elongating cases. Still, this effect would have close to zero bias over long datasets, when the Wave Glider navigates both opposing and following the waves.
9) Determination of bulk wave parameters
Bulk wave parameters are calculated from the spectral results, including significant wave height , peak wave period , energy period , and peak wave direction . As the ensemble spectra are calculated every 30 min, these bulk parameters are also determined every 30 min. Wave heights cannot be directly calculated from the raw data without first computing the spectra, because the conversion from horizontal velocities to wave amplitude requires knowledge of the wave frequencies [see Eq. (9)]. The significant wave height, also often referred to as , is defined as
The peak wave period is the inverse of the frequency at the maximum of and the peak wave direction is the value of at the maximum of . The exact peak of the spectral energy is often poorly defined, given that typical open ocean conditions have multiple wave systems present (i.e., mixed seas). Thus, a weighted energy period, defined as
is often preferred. In the present work, both estimates of the bulk wave period are evaluated.
10) Evaluation datasets (wave)
Several datasets have been collected to evaluate and characterize the wave measurements and the onboard processing routines. There are five datasets:
Southern Ocean (Drake Passage), winter 2017
Astoria Canyon (offshore Oregon), winter 2017
Hawaii, autumn 2016
Monterey Canyon (offshore California), spring 2015
Hawaii, spring 2014
For each of the datasets listed, groundtruth data are provided by an independent sensor, either a Datawell MOSE-G unit on board the vehicle or a Datawell Waverider buoy (part of the CDIP network) within 10 km of the vehicle. All datasets are included in the overall evaluation given in section 3. The Southern Ocean and Astoria Canyon cases are given extra attention, as these are the most extensive datasets. The more recent of the Hawaii datasets is also of particular interest, as this evaluates the new Microstrain GX4 sensor as the source of GPS and IMU data (the actual variables of are the same as from the GX3). Figure 2 shows the time series of bulk wave parameters from the onboard GPSwaves results and the Datawell reference for both the Southern Ocean and Astoria Canyon datasets. The Astoria Canyon measurements are the preferred evaluation data, because the ground truth data are fully independent (CDIP station 179) of the vehicle motion. However, the Southern Ocean dataset is significantly longer in duration and has a full range of wave directions (whereas the Astoria Canyon directions cover only 40°–330°). The ground truth data for the Southern Ocean dataset are provided by a Datawell unit on board the vehicle.
b. Wind measurements and processing
The approach to estimating wind stress from the Wave Glider follows closely the work of Edson et al. (1991) and Yelland et al. (1994), in which fast-sampling three-axis sonic anemometers are used to measure turbulent fluctuations of the wind velocity components in the inertial subrange of the turbulent length scales. The essential point is that the inertial subrange is observed at frequencies beyond the influence of platform motion. After using this portion of the frequency spectrum to estimate the turbulent dissipation rate ε, the wind friction velocity is inferred by assuming a “law of the wall” balance of production and dissipation. The wind stress is then , where kg m−3 is the density of air. This method captures only the turbulent stress; there is evidence that additional wave-supported stresses are important in the atmospheric boundary layer (e.g., Hara and Sullivan 2015) but that contribution to the wind stress is not addressed with these methods.
1) Raw wind data collection
Raw turbulent winds are collected with a Gill Instruments WindMaster three-axis ultrasonic anemometer (Hampshire, United Kingdom) mounted to the bow of the Wave Glider at 0.77 m above the sea surface. The Gill samples at 10 Hz continuously; the raw data of relative wind velocity components are stored on board the vehicle and processed every 30 min. (note the similarity to raw wave data in both notation and burst interval). The Gill is oriented such that the negative u velocity component is a relative wind coming at the bow of the vehicle. Reference wind data are collected using an Airmar WX200 weather station (Milford, New Hampshire) mounted in the center of the vehicle at 1.05 m above the sea surface. The Airmar collects two-axis (horizontal) wind data at 1 Hz and performs onboard corrections for pitch, roll, and heading to estimate the true wind; the raw data are not retained and a 10-min average value for true wind speed and direction is reported.
2) Adjustment to 10-m wind speed
Both measurements of the wind are adjusted to a 10-m reference height, assuming a logarithmic profile (e.g., Zedler et al. 2002). The Gill measurements of relative wind are first adjusted for burst-averaged vehicle heading and speed over ground to produce burst-averaged true winds. Raw motion correction is not attempted, because the raw data recording is not synchronized with the Microstrain GPS/IMU data recording. At moderate wind speeds, this height adjustment is approximately for the Gill and for the Airmar, where the capital indicates the burst-average wind speed (i.e., the magnitude of the horizontal components ).
3) Turbulent wind spectra
Calculation of the turbulent wind frequency spectra follows closely the steps for spectral processing of the wave data but with Hz and s. An additional step for the wind spectra is the merging of every 11 neighboring frequency bands before ensemble averaging the windowed spectra. This step increases the degrees of freedom and smoothness of the resulting ensemble spectra, at the expense of a lower-frequency resolution. Spectral densities are estimated for each component of wind velocity, and these are expected to have a shape in the inertial subrange, which is a manifestation of the wavenumber energy cascade (Taylor 1937; Kolmogorov 1941). The inertial subrange is typically in the frequency range of 1–5 Hz, although it may be obscured by instrument noise at high frequencies (Edson et al. 1991). Here instrument noise is assumed to be white (i.e., uniform in frequency), and an a priori estimate of a m2 s−2 Hz−1 noise floor is subtracted from all frequencies of each spectrum.
4) Quality control of wind spectra
The turbulent wind spectra are used only when the Wave Glider is head to wind, such that the wind measurements are not in the wake of the other masts and antennas on the vehicle. This is determined by screening the burst-average components for cases when , which indicates the relative winds are within of the bow. Navigational constraints and other mission objectives may make this difficult to achieve more than a fraction of the time, but it is essential for the quality of the results. Next, the wind spectra are screened for isotropy in the inertial subrange, using the ratio of the spectral levels of the horizontal components to the spectral levels of the vertical component. This parameter would be equal to 1.0 for perfect isotropy; a maximum anisotropy of 1.5 is allowed. Finally, the wind spectra are screened for the expected shape in the inertial subrange, using the ratio of the standard deviation to the mean of the magnitude of the compensated spectrum . This quality parameter is equal to 0.0 for a perfect theoretical spectrum; a maximum quality ratio of 0.6 is allowed.
5) Wind stress from inertial dissipation
The turbulent wind spectra are averaged as compensated spectra over the range Hz to determine the dissipation rate,
where is the Kolmogorov constant and is the burst-averaged advection speed prior to the adjustment to the reference height of 10 m. The frequency range Hz is selected as a portion of the spectrum that is above wave motions and wave propulsion contamination, yet below significant instrument noise. Assuming neutral atmospheric stability and a constant stress layer, the wind friction velocity is
where is the von Kármán constant and m is the measurement height. The wind stress is simply
which can be related to standard drag laws using
where is the nondimensional drag coefficient.
6) Evaluation datasets (wind)
The wind stress method was only recently implemented on the Wave Glider, so the Southern Ocean dataset is the only case available for evaluation. The dataset is 68 days in duration, and the burst-averaged true wind speeds range from 1 to 17 m s−1. Time series of wind magnitude and direction are shown in Fig. 3. The screening for head-to-wind conditions is severe in this dataset; only 30% of the 30-min bursts meet the criterion. These are mostly from the first half of the mission. During the second half of the mission, the vehicle was on a downwind crossing of the Drake Passage, and the winds were almost aways from astern. Note that this also corresponds to the poor estimates of wind direction from the Gill measurements in Fig. 3b.
The Southern Ocean dataset lacks a fully independent wind measurement, such as a nearby moored buoy, to evaluate the application of the inertial dissipation method on the Gill data. There are 2 h of validation data available from a test mission on the Washington coast, in which results from a nearby direct covariance wind system (Edson et al. 1998) on the R/V Jack Robertson at 6-m height will be compared with the Wave Glider methods. Otherwise, reference data will use the bulk drag formula applied to the Airmar wind measurement and a wind reanalysis product from the U.S. National Centers for Environmental Prediction (NCEP; Saha et al. 2010). Finally, the wind-wave equilibrium stress [Eq. (1)] will be used as another reference and as a demonstration of the combined wave and wind stress results.
a. Wave results
1) Wave spectra
Since the determination of wave spectra is a precursor to the determination of bulk wave parameters in the method, the wave spectra results are also presented first. Average spectra are shown for the Southern Ocean and Astoria Canyon datasets. These averages are not valid wave spectra, in the sense of having stationary statistics over such a wide range of wave conditions; however, they are illustrative of the signal characteristics and overall quality. Averages of the ratios between energy spectra are also presented; these are more robust, as they are relative metrics that are calculated for each 30-min burst of data (and then averaged over the whole dataset).
2) Spectral wave energy
The averaged wave energy spectra in Fig. 4 show good agreement between the GPSwaves method and the Datawell reference. There is some bias in the range of approximately Hz, where GPSwaves estimates more energy than the Datawell. This is most evident in the average ratio of the energy spectra, , which is computed for each pair of spectra and then averaged over the whole dataset. If the directional spread screening of was not done, then there would also be a bias at low frequencies. The bias is broadband, and thus it does not appear related to a particular wavelength, as might be expected if there were contaminations related specifically to the tether length of the Wave Glider (the tether length was 8 m for the Southern Ocean mission and 16 m for the Astoria Canyon mission).
The bias in spectral energy is likely related to wave propulsion, which would contribute nonorbital motion to the raw velocity signals. Removal of this signal would require detailed knowledge of the vehicle motion to develop a transfer function that could correct for these contaminations. Alternatively, an empirical spectral correction factor could be defined based on the observed spectral ratios (i.e., Figs. 4a,b). Although it is tempting to apply such a correction, , the empirical nature could cause significant errors when applied to conditions outside of the training datasets. This is especially risky for the case of small waves. As shown in the next section, the net effect of the contamination is less than a 5% bias high in , and thus for many applications it may be prudent to use as reported.
3) Spectral check factors
A commonly used quality metric for wave energy spectra is the “check factor,” which is the ratio of vertical motion to horizontal motion at each frequency band (e.g., Thomson et al. 2015). In deep water, this value is theoretically equal to one, because wave orbits are circular. Figure 5 shows the average check factors from the Astoria Canyon dataset. Reference check factors are not shown for the Southern Ocean dataset, because they are produced only for the CDIP Datawell buoys and not from the Datawell units on board the vehicle (so there are no reference check factors for the comparison in the Southern Ocean dataset). As expected, the Datawell reference check factors are close to one for the frequencies dominated by wave energy, and deviations from one are seen only at the very lowest or highest frequencies. The GPSwaves results, by contrast, are less than one at most wave frequencies. This indicates that horizontal motions exceed vertical motions, consistent with the mechanism listed above for a bias high in spectral energy: wave propulsion. Again, it is tempting to use the spectra check factor to correct for the bias, much as with the that could be determined by the spectral energy ratio, but such an empirical approach built from average conditions may be ineffective (or worse) when applied to other conditions.
4) Spectral wave directions
The averaged spectral wave directions in Fig. 6 show excellent agreement between the GPSwaves method and the Datawell reference. The only deviations are at the lowest frequencies, where the wave signals are generally very small. The average ratio between the individual estimates is close to one, except at the higher frequencies. This suggests that the directional wave estimates from Wave Glider may actually be better than the scalar energy estimates—at least there is no obvious bias. This is likely because the directional moments are normalized quantities (see the appendix), and the fidelity of the estimates comes from the quality of the relative phases of the raw signals.
The directional spread also agrees well across all frequencies (not shown). When compared with the independently moored Datawell, as in the Astoria Canyon dataset, the GPSwaves results have spreads that are biased high by 5°–10°. This bias is not observed when comparing with the Datawell unit on board the vehicle, as in the Southern Ocean dataset, and the implication is that vehicle motion broadens the directional estimates slightly.
An example of the directional moments () are shown in Fig. 7 for the randomly selected case at 0000 UTC 29 January 2017 from the Astoria Canyon mission. The corresponding directional spectra are estimated following the MEM of Lygre and Krogstad (1986) and shown in Fig. 8. The GPSwaves results from on board the vehicle are in good agreement with the Datawell Waverider moored at CDIP station 172. There are some subtle differences at higher frequencies, but the larger pattern is well captured. A comparison of all spectral moments indicates the relative errors are typically 5% or less with a bias of less than 2%.
5) Bulk wave parameters
The bulk wave parameters , , and are compared as scatterplots in Fig. 9 for both the Southern Ocean and the Astoria Canyon datasets. The bulk wave parameters , , and are compared for all evaluation datasets in Tables 1–3. Two skill metrics are calculated: the bias and the root-mean-square (rms) error. The bias is signed using the convention of the GPSwaves result minus the Datawell reference, so a positive bias indicates that the GPSwaves result tends to be higher than the Datawell. The rms error is symmetric (i.e., it is calculated after determining and removing the bias), and it indicates the level of uncertainly in the GPSwaves result. It should be noted that some reference data are unavailable, because the Datawell units had dropouts and sample gaps (368 bad bursts out of 3096 total, or 11%, for the Southern Ocean dataset).
6) Wave height comparison
The significant wave heights are generally in agreement between the onboard GPSwaves result and the Datawell reference. This is expected from the wave energy spectra comparisons, since the wave height is simply 4 times the square root of the integrated spectra [Eq. (12)]. Both the bias and the error are larger for the datasets with larger waves. The relative skills are, conservatively (worst case), 2% bias and 5% symmetric error. For full ocean conditions, this indicates that the GPSwaves method running on board the Wave Glider provides wave height estimates to within cm (about the length of a standard desk ruler).
7) Wave period comparison
There is a strong negative bias in the estimates of peak wave period from the GPSwaves method on the Wave Glider. This is a direct result of the bias in the wave energy spectra, which has more energy at higher frequencies than the Datawell reference spectra. The negative bias is up to 25% of the range, while the symmetric error is up to 30%. These errors are exacerbated by the noisy nature of the peak period metric itself and are further enhanced by the inverse relation of period to frequency (i.e., the uniform frequency bands from the FFT give nonuniform spacing of the wave period with especially poor resolution at low frequencies).
Many researchers prefer an energy-weighted wave period [Eq. (13)] because it is a more stable estimator. These values compare much more favorably between the GPSwaves results and the Datawell reference. For the case of the Southern Ocean, the wave period bias reduces from −1.5 to 0.1 s, and the rms error reduces from 2.1 to 0.8 s. Thus, a more stable measure of the wave period can be reproduced by GPSwaves with less than 1% bias and 5% symmetric error.
8) Wave direction comparison
Peak wave directions agree well between the GPSwaves method on the Wave Glider and the Datawell reference with small bias. There is, however, large scatter, reflected in rms errors up to 34°. Much of this scatter is derived from disagreements in the peak wave period, which sets the frequency at which the dominant direction is selected from . If only results with good agreement in are used, or if the value of is selected consistently based on from both data products, then much better agreement is found. The errors for these cases are reduced to less than 15°, and the bias remains small. The corresponding relative skill metrics are 1% negative bias and 4% symmetric error.
b. Wind results
1) Bulk wind speed and direction
The bulk wind speed and direction agree well, although the direction from the Gill is error prone if the Wave Glider is not pointed into the wind (see right half of Fig. 3b). As the Airmar is on the tallest mast, above all obstructions, it does not have this limitation. The adjusted 10-m wind speeds from the two instruments are directly compared in Fig. 10, with a bias of −0.2 m s−1 and a symmetric error of 1.3 m s−1 reported in Table 4. The wind speeds are a relative bias of 1% and relative error of 8%. This agreement in from instruments at two different heights suggests validity in the log-layer adjustment. This is further validated by the agreement with the NCEP reanalysis winds (also at 10-m reference height). For the higher wind speeds, the Gill measurements are somewhat low relative to the Airmar (a bias of −0.8 m s−1 and an error of 1.4 m s−1 for times when m s−1), and this may be the result of wave sheltering.
2) Turbulent wind spectra
The turbulent wind spectra are shown in Fig. 11, using wind speed bin averages. As expected, the spectral levels are sorted according to the bulk wind speed. All wind spectra show a motion-contaminated region from to Hz, which is consistent with the peaks in the wave energy spectra (Fig. 4a). Above these frequencies, there is an inertial subrange following . This range is used for estimation of the dissipation rate [Eq. (14)] and the subsequent estimation of the wind friction velocity [Eq. (15)]. These spectra are corrected for white noise, which would otherwise flatten the spectra at the very highest frequencies.
3) Wind stress
The wind friction velocity estimates are shown in Fig. 12 and Table 4, with a comparison to a standard drag law (Smith 1980) and to the wave equilibrium estimates [Eq. (1)]. The wind stress is simply the square of these estimates [Eq. (16)]. There is good agreement with 0.01 m s−1 bias (1%) and 0.03 m s−1 symmetric error (4%) relative to the drag law. There is also good agreement for two additional validation points; these compare with the direct covariance system on the R/V Jack Robertson during a test mission on the Washington coast in July 2016. There is increased variability in the results at higher winds, which may be related to changes in wave age and or changes in surface roughness as a result of wave–current interactions in the Drake Passage.
A stability correction can be made between Eqs. (14) and (15), but the effect is typically only a few percent change in . For the Southern Ocean dataset, air temperatures were always warmer than ocean temperatures, by 2°–5°C, and the correction for these stable conditions was less than 1% throughout.
a. Lack of propulsion contamination
As noted, the wave height bias is likely related to wave propulsion, which raises the question: How can a Wave Glider measure waves to within 5% while simultaneously using the waves for propulsion? The heuristic answer is in the relative scales of motion and the inefficiency of the propulsion. First, the Wave Glider advances forward only a few meters with the passage of each wave, such that the raw signal of the horizontal motion (and the variance) is dominated by the natural wave signal. Second, the surface float is towed by the sub, not the other way around (i.e., the sub is the part that is restricted from moving backward after a wave passes and thus gives propulsion). This means that the surface float is free to track with the wave orbitals at the surface with the weak constraint of being pulled along by the sub. This is not all that different from a mooring tether keeping a CDIP buoy on station in a strong ambient current (although CDIP moorings often employ a rubber cord to soften this constraint). The same is true for moored air–sea flux buoys that estimate wind stress while moving with passing waves.
b. Natural wave variability
For the datasets at Astoria and Monterey Canyons, which have moored CDIP buoys as the Datawell reference, it is important to note the known natural spatial variations in wave fields. This would be a secondary explanation for the difference in wave results. As shown in Gemmrich et al. (2016), the separation of up to 10 km between the CDIP mooring and the Wave Glider can result in a 5% difference in wave heights. Similarly, there is significant uncertainty in determining bulk wave parameters from finite-length records. Drazen et al. (2016) show that the degrees of freedom typical for processing 30-min bursts of wave data can result in significant wave height uncertainties of up to 20%. Although these combined effects can account for the errors in the comparisons, they are less likely to account for the biases.
c. Alternative approaches to wind stress estimation
Following the work of Edson et al. (1998), many researchers now use a direct covariance method to calculate the wind friction velocity as , where are the turbulent fluctuations in the wind velocity components, after removing platform motion and mean winds. This was the method used for the two points of a ship comparison in Fig. 12. For the Wave Glider method, motion is not measured synchronously, nor collocated with turbulent winds, and thus motion correction has not been successful. The inertial dissipation method used herein is far more robust to motion contamination, because the method isolates very small scales in the turbulent wind field. Another option might be to fit Kaimal et al. (1972) curves to the cospectra of at higher frequencies and allow the fitting to span across the lower and wave-contaminated frequencies but that has not be applied here.
Datasets spanning full open ocean conditions demonstrate that measurements of wave directional spectra and wind stress, along with bulk wave parameters and bulk winds, can be accurately made from the Wave Glider platform. In an ongoing upgrade, raw GPS and IMU data will be collected with a Microstrain GX4 sensor, rather than a Microstrain GX3 sensor, and a subset of the data confirms that a new sensor will supply raw data of similar quality for the GPSwaves algorithm. Future work on this topic would develop a spectral transfer function for the Wave Glider, such that the platform response and propulsion could be directly removed from wave and wind spectra.
Thanks to APL field engineers Alex de Klerk, Avery Snyder, and Joe Talbert for their help with the preparation, testing, and deployment of the Southern Ocean Wave Glider. Thanks to the captains and crews of all ships involved in these tests. Thanks to CDIP for providing excellent ground truth data from a consistently maintained database with publicly available software and documentation. Thanks to Liquid Robotics for integration and general support. Thanks to Seth Zippel and two anonymous reviewers for comments that improved the manuscript. Funding for the Southern Ocean work was provided by the National Science Foundation (PLR1558448). Funding for the other datasets provided by Liquid Robotics, Inc. (producers of the Wave Glider). Data are available online (http://faculty.washington.edu/jmt3rd/Waveglider/).
Directional Moments from Cross Spectra
Directional moments from the cross-spectra are defined in a Cartesian system relative to east, consistent with raw data that are velocity components of u (positive east) and υ (positive north). The first two moments use the cross-spectra of the velocity components and the vertical acceleration,
and are normalized by the corresponding autospectra such that they vary from −1 to 1. The horizontal velocity and vertical acceleration of propagating waves are in phase, and thus the real part, or “cospectrum,” of the signals is used. This is an important distinction from the conventional definition, which uses signals from horizontal and vertical displacements or from surface slopes and vertical displacements. Those signals are out of phase, and thus the imaginary part, or “quad spectrum,” of the signals is used in the conventional definitions (Kuik et al. 1988; Herbers et al. 2012).
Figure A1 shows how the values are used to identify the principal axis of wave motion at each frequency. The dominant wave direction and spread at every frequency are given by
These are the directions the waves are coming from a clockwise value relative to east. This must be converted to counterclockwise (change sign) and made relative to north (add 90°) to get a direction in the nautical convention. This direction is already true, not magnetic, as the raw data are already relative to true north (and true east). This is another advantage of using data that are not dependent on heading data (and therefore not dependent on the quality of magnetometer calibrations).
The other two moments use the cross-spectra of the velocity components with each other,
and are again normalized by the corresponding autospectra such that they vary from −1 to 1. The horizontal velocities are in phase with themselves and each other, and thus the autospectra and the real part of the cross-spectra are used.
Figure A2 shows the how the values are used to identify the principal axis of wave motion at each frequency. The dominant wave direction and spread at every frequency are given by
These directions are again clockwise values relative to east, but they cover only two (out of four) quadrants. There is a 180° ambiguity in the result, such that can distinguish only the principal axis along which the waves are propagating but not whether waves are to or from. Thus, is the preferred estimate of wave direction, though can be used when near a coastline (where waves must be propagating from offshore, and the ambiguity can be resolved by common sense). The term is also useful when vertical acceleration data are not available, because depend on the velocities alone.