1 Introduction

Saving lives at sea and on beaches is a United States Integrated Ocean Observing System (IOOS) priority that is supported by a Memorandum of Understanding between the National Oceanic and Atmospheric Administration (NOAA) and the United States Coast Guard (USCG). Nationally, the Coast Guard saves an average of 14 lives each day. Unfortunately, another three lives a day are lost (http://www.uscg.mil/hq/cg5/cg534/sarfactsinfo/USCG_SAR_Stats.asp). To reduce the lives lost, the critical USCG need is to optimize search and rescue (SAR) operations to minimize search time (Arthur Allen, personal communication). The USCG decision support tool for search and rescue cases is the Search and Rescue Optimal Planning System (SAROPS). SAROPS uses observed or predicted surface wind and surface current fields from the USCG’s Environmental Data Server (EDS) to predict the trajectories of floating objects. During an actual event, a cluster of typically 5,000 virtual objects is deployed in the EDS provided surface wind and current fields and allowed to drift over time (Spaulding et al. 2006). Based on the distribution of the predicted cluster, the tool helps search planners coordinate the appropriate search strategy given available assets. Environmental input fields include in situ and remote observations and numerical and statistical forecasts.

The USCG and the ocean observing community have worked to evaluate the value of high frequency (HF) radar-derived surface currents within SAROPS. These studies have shown that the time series of mapped surface currents provided by regional scale HF radar networks reduce search areas by a factor of 3 (O'Donnell et al. 2005; Ullman et al. 2003; Roarty et al. 2010). If these surface currents have lower uncertainties, there is less dispersion in the cluster, a smaller more accurate search area, and greater likelihood for success.

Coordinated through the U.S. IOOS, the Mid-Atlantic Regional Association Coastal Ocean Observing System (MARACOOS) has evaluated the application of these data and provided estimates of uncertainty for incorporation into SAROPS (O'Donnell et al. 2005; Roarty et al. 2010; Ullman et al. 2003). The 35-site MARACOOS network is a triple-nested HF radar network that covers the continental shelf between Cape Hatteras, NC and Cape Cod, MA with higher resolution coverage nearshore and into the estuaries. Since the first sites were deployed in 1998, there have been evaluations to determine the optimal system configuration and associated vector uncertainty in the 5-, 13-, and 25-MHz bands (Kohut and Glenn 2003; Kohut et al. 2006; O'Donnell et al. 2005; Roarty et al. 2010). These frequencies cover the scales of the mid and outer shelf, inner shelf, and harbors/estuaries, respectively. Based on this work, in May of 2009, the MARACOOS HF radar network became the first to support the operational SAROPS decision support tool. To maintain this operational status within SAROPS, the surface current data and associated forecasts must continue to meet quality standards defined by the Coast Guard.

In the USA, there is now an effort to expand the operational support of SAROPS beyond the Mid-Atlantic Bight to the national scale (Harlan et al. 2010). In that process, there is currently an evaluation to determine the optimal processing suite to maintain a national surface current mapping product that meets USCG quality standards. In this context, we use the HF radar network in the Mid-Atlantic Bight (MAB) to evaluate two different algorithms that both combine radial component vectors from an array of coastal stations into hourly total vector surface current maps. This evaluation is enabled through a partnership between the USCG Research and Development Center, the USCG Office of Search and Rescue, and MARACOOS. It is based on two available combination algorithms, the established unweighted least squares (UWLS) (Lipa and Barrick 1983) already providing surface currents to the EDS to support operational SAR in the MAB and the more recently adapted optimal interpolation (OI) (Bretherton et al. 1976; Kim et al. 2007). There are two main objectives of this study. The first is to evaluate the recently adapted OI-derived surface current maps with parameter inputs consistent with the existing UWLS implementation. The control for this evaluation will be the already implemented UWLS product. Since the OI implementation is new to the HF radar community, the second objective is to quantify the sensitivities of the OI-derived surface current maps to varying parameter inputs. For this objective, the base OI solution evaluated in objective 1 will serve as the control and be compared to other OI solutions with varying spatial decorrelation scales and error thresholds. Both objectives will be discussed in the context of USCG search and rescue using in situ current estimates and drifter trajectories. Background information and specific methods used are discussed in Sections 2 and 3. Expected differences due to known environmental variability in the surface current fields and evaluation results are presented in Sections 4 and 5. Finally, Sections 6 and 7 discuss the results and general conclusions.

2 HF radar processing

2.1 HF radar surface current processing

HF radar systems, typically deployed along the coast, use Bragg peaks within a transmitted signal (3~42 MHz) scattered off the ocean surface to calculate radial components of the total surface velocity at a given location (Barrick et al. 1977). Peaks in the backscattered signal are the result of an amplification of a transmitted wave, at grazing incidence, by surface gravity waves with a wavelength equal to half that of the transmitted signal (Crombie 1955). The frequency of the backscattered signal will be Doppler-shifted depending on the velocity of the scattering surface. Using the linear wave theory, the phase speed of the surface waves can be separated from the total frequency shift, leaving only that shift due to the surface current component in the direction of the antenna. Over a given time period, sites along the coast generate radial maps of these component vectors with typical resolutions on the order of 1–6 km in range and 5° in azimuth. The HF radar sites in the MARACOOS network are all SeaSonde direction finding systems manufactured by Codar Ocean Sensors (Barrick 2008). The direction finding radars use a three element receive antenna mounted on a single post to determine the direction of the incoming signals. The angular resolution, set in the processing, is 5° (Barrick and Lipa 1996; Teague et al. 1997). Since the Doppler shift can only resolve the component of the current moving toward or away from the site, information from at least two sites must be geometrically combined to generate total surface current maps.

The main source of error in the geometric combination from radials to totals is based on the geometry of the network (Chapman et al. 1997). The contribution of the geometric error is based on the relative angles between the radial component vectors at the total grid point, commonly referred to as the geometric dilution of precision (GDOP). The geometric error increases as the angle between radials from contributing sites moves away from orthogonality.

In HF radar networks operating around the world for the past few decades, the most commonly used algorithm to combine maps of radial velocities into total vector currents is UWLS (Lipa and Barrick 1983). The coverage of the UWLS surface current maps is dependent on the density of available radial components measured at each site. Due to variations in signal-to-noise based on either the electromagnetic interference or sea state conditions, the available radial data can vary in both space and time. To overcome these data gaps, various interpolation techniques have been applied to HF radar total vector fields. These algorithms that include 2DVAR (Yaremchuk and Sentchev 2009), normal modes (Lipphardt et al. 2000), open modal analysis (Kaplan and Lekien 2007), and statistical mapping (Barrick et al. 2012; O’Donnell et al. 2005) have largely been applied to the UWLS surface current maps after the UWLS combination. Recently, Kim et al. (2007) introduced a method that interpolates data as part of the combination step from radial component vectors to total vector maps. This approach weights radials based on the decorrelation scale of the observed current fields. Details on both the UWLS and OI algorithms as implemented in this study are described in Sections 2.2 and 2.3.

2.2 Unweighted least squares

UWLS is the most commonly used algorithm to combine fields of radial surface current components from a network of shore sites into a total vector map over a fixed grid (Lipa and Barrick 1983). The UWLS calculation provides a total current vector that minimizes the error between the radial components of the calculated total vector and the measured radial velocities. The algorithm assumes that all vectors within a defined area relative to the grid point have an equal weight in the calculation. The UWLS makes no assumption on the signal variance of the surface currents and assumes that the correlation function of the surface currents within this radius is unity and zero outside. The contribution to the GDOP error in the implementation of the UWLS is quantified as a mapping error (Gurgel 1994). For each total vector estimate, there is an associated mapping error value that is scaled based on the contribution to GDOP. The larger the mapping error value, the more significant the error due to geometry is in the resulting total vector. For this study, and consistent with the UWLS implementation in the national network for the 5-MHz radial data, total vectors were calculated with radial components that were within a radius of 10 km from each grid point and a mapping error threshold of 1.5. If any UWLS-derived total vector had an error threshold above 1.5, it was not considered in the analysis. This is consistent with quality standards already governing the data going to the EDS and SAROPS and the national network (Harlan et al. 2010; Roarty et al. 2010).

2.3 Optimal interpolation

The technique for the optimal interpolation of oceanographic data was derived from the Gauss–Markov theorem (Bretherton et al. 1976). This technique was more recently adapted to calculate total vector surface currents from available HF radar radial component observations (Kim et al. 2007). The OI requires a quantification of the surface current signal variance and error measurement variance. For our study, the signal variance is the expected variance of the observed surface current field and the error measurement variance is the uncertainty assigned to the input radial velocities. It also requires a decorrelation length scale in the east and north direction to define the correlation function of the surface currents. We define a base OI run for this study that uses a defined decorrelation scale of 10 km. This is consistent with the scale of the radius used as a community standard in the UWLS algorithm.

Radial velocities are weighted based on an exponential decay defined by the decorrelation scales in the north (S y ) and east (S x ) directions. Vectors closer to the grid point are weighted more than those further away. S x and S y set the length scale of the decay in this weight. Radials that fall outside this scale are weighted much less than those closer to the grid point. Unlike the UWLS in which all radials within the radius are weighted 1 and included in the calculation, the OI can use vectors outside the defined length scale. For example, if a radial vector were 10 km from the grid point, it would be weighted as 1 for the UWLS and 0.37 for the OI. Another vector 20 km from the grid point would be weighted 0 for the UWLS and 0.14 for the OI. Therefore, the OI will use radial data from further away than the UWLS, but these vectors will be weighted significantly less. The other input parameters for the OI algorithm are (1) the signal variance of the surface current fields and (2) the data error variance of the input radial velocities. The signal variance for all OI solutions was set to 420 cm2/s2 for all grid points in the domain; 420 cm2/s2 is representative of the typical conditions in the MAB with surface currents on the order of 20 cm/s (Beardsley and Boicourt 1981; Dzwonkowski et al. 2009, 2010; Gong et al. 2010; Kohut et al. 2004). The data error variance is defined by Kim et al. (2007) as the sum of the average measurement uncertainty and the average standard error of the radial currents (Eq. 17, Kim et al. 2007). These two values are calculated from the input radials. The average measurement uncertainty is the mean of the uncertainty values reported with each radial velocity and the average standard error is calculated from radial velocities themselves. Applying Eq. 17 from Kim et al. (2007) to the radial data and associated uncertainties from each of the four sites considered in this study, the data error variance is 66 cm2/s2 and held constant across the field. All OI solution sets evaluated in this study use a constant 420 and 66 cm2/s2 for the signal variance and data error variance, respectively.

To test the sensitivity of the OI solution to varying decorrelation scales, we considered a second parameter set that stretched S y along the isobaths. In the MAB, the decorrelation scales tend to be longer along the isobaths (Beardsley and Boicourt 1981; Kohut et al. 2004). Therefore, the decorrelation scale for the second run was elongated by a factor of 2.5 and rotated along the isobaths (S x  = 10; S y  = 25). To do this, we rotated the coordinate system that defines S x and S y in the along-isobath direction and assigned the along-isobath component (S y ) a value of 2.5 times larger than S x . The algorithm uses these parameters to determine the weight of the individual radials going into the combination. So the result is a weight based on decorrelation scales oriented along isobaths. The along-isobath direction was defined as 31° clockwise from true north, consistent with previous generalizations of the topography of the central MAB (Kohut et al. 2004).

With the OI algorithm, the geometric contribution to the error of the total vector estimate is quantified as individual error estimates for each component. For each component of the total current vector, there is an associated normalized uncertainty. The threshold for the base OI run was set to 0.95. This threshold was selected to be consistent with the implementation of the OI algorithm in the Southern California Bight (Sung Yong Kim, personal communication). The Southern California system was the only known real-time implementation of the OI algorithm at the time of our study. If any total vector had an uncertainty greater than 0.95 in either the east or north component, it was not considered in the evaluation. The sensitivity of the OI results to this threshold was tested with a lower value of 0.6. This second value was selected based on the balance between data quality and data coverage (Section 5.2 of this manuscript). In summary, there will be four OI runs tested. The first with a decorrelation defined as S x  = S y  = 10 and an error threshold of 0.95 will be compared to the UWLS solution with a radius of 10 km and a mapping error threshold of 1.5 to meet our stated objective 1. For objective 2, an additional three OI runs with varying decorrelation scales (S x  = 10; S y  = 10; and S x  = 10, S y  = 25) and error thresholds (0.95 and 0.6) will quantify the sensitivity of this recently implemented combination algorithm to these settings.

3 Data

3.1 HF radar

The HF radar data were provided by the MARACOOS network, a nested array of 5 MHz long range systems, 13 MHz standard range systems, and 25 MHz high resolution systems (Roarty et al. 2010). The nested coverage stretches across the shelf from Cape Cod, MA to Cape Hatteras, NC with higher resolution into the estuaries. This study focuses on the radial data from four 5 MHz sites in the central MAB in Sandy Hook, NJ; Loveladies, NJ; Wildwood, NJ; and Assateague, MD between February 20, 2007 and April 30, 2007 (Fig. 1). Each remote site was operated with the Quality Assurance/Quality Control recommendations from the MARACOOS operators and the Radiowave Operators Working Group community. These are the same data provided to the national HF radar server at the NOAA National Data Buoy Center (http://hfradar.ndbc.noaa.gov/). Every hour, the available radial velocities are combined into a single total vector map on the national network 6-km grid (Temll et al. 2006). A total vector was only generated if at least three radial velocities from at least two remote sites were available to the combination algorithm. For each grid point, the available radial velocities within the defined area were combined into a single total vector with both the UWLS and OI algorithms described in Section 2. For the comparison to the acoustic Doppler current profiler (ADCP) moorings, the UWLS and OI totals were calculated on a grid that exactly matched the mooring locations. For drifter estimated current and trajectory comparisons, the OI and UWLS totals calculated on the national 6-km grid were spatially interpolated to each drifter position at each time step.

Fig. 1
figure 1

Map of the study area showing the location of the HF radar stations (circles), the ADCPs (stars), and track of the February–March drifter (thick black line) and April drifters (dashed lines). The isobaths (meters) are shown as faint gray lines. All drifter tracks start at the triangles and end at the squares. The total vector coverage for the UWLS solutions is shown for the 50 % (thick black) and 90 % (thin black) contours

Full size image

3.2 Acoustic Doppler current profilers

Four ADCPs were deployed off the coast of New Jersey as part of the National Science Foundation supported Mid-Shelf Front Experiment (Ullman et al. 2012). They were deployed on January 13, 2007 and recovered 83 days later on April 6, 2007. Three of the moorings were deployed in a cross-shelf line approximately 10 km apart (Fig. 1). A fourth mooring was deployed 11 km upshelf of this cross-shelf array. Since the offshore mooring was recovered from a different location than the original deployment, special processing steps were taken with the current profile time series. The data record for this mooring indicated that on February 12th, about 32 days into the deployment, there was an hour of bad data. This was followed by 53 days of good data. Scrapes on the unit seen after recovery indicate a likely strike, staying upward looking throughout the event. Based on the deployment and recovery locations, it is estimated that the ADCP moved approximately 2.4 km. Since this is within the approximate 20-km scale of the HF radar measurement cell, these data were merged and treated as one time series. The position was the duration-weighted mean of the deployment and recovery locations. All units were configured with 1-m bins in the vertical. Sampling was set to collect a 10-min ensemble each hour (Codiga 2007). These hourly data were averaged into 3-h files outputted every hour to be consistent with the sampling and averaging of the HF radar processing.

3.3 Drifters

All drifters used in this evaluation were self-locating data marker buoys (SLDMBs) provided by the USCG Research and Development Center. These drifters are drogued to 1-m depth with minimal surface expression to minimize wind drift (Allen 1996). The first deployment included two drifters released on February 24, 2007. Of the two drifters deployed, only drifter 43484, deployed inshore of the Mid-Shelf front (Ullman and Cornillon 1999), remained in the HF radar coverage (Fig. 1). This winter deployment overlaps the ADCP deployments providing a contrasting Eulerian and Lagrangian evaluation in regions of the HF radar footprint beyond the reach of the ADCP array. The second deployment of six drifters began April 3, 2007. These six drifters remained primarily in the HF radar data footprint through the end of April (Fig. 1). While this second group of drifters did not overlap with the ADCP deployment, the larger cluster of drifters provided simultaneous observations across different regions of the HF radar data footprint. Throughout the deployment, velocities based on two drifter positions 1 h apart were calculated every half hour. Drifter-derived currents estimated every 30 min were averaged into 3-h files outputted every hour to be consistent with the sampling and averaging of the HF radar processing. Three of the drifters entered the Gulfstream toward the end of the time series (43266, 43271, and 43325). For these three drifters, only the position data over the shelf were included in the evaluation.

3.4 Drifter trajectories

Drifter tracks were used to evaluate particle trajectories estimated from each algorithm. Each hour, we took the actual position of drifter #43484 as the initial condition for the particle release. This drifter was selected because most of its path is within the 90 % coverage contour of the HF radar data giving a longer time period with which to evaluate the trajectories (Fig. 1). The entire track of the drifter was split into 713 forty-eight-hour segments each beginning 1 h after the last. For each segment, we released 1,000 particles into the HF radar fields at the location of the actual drifter at time 0. The particles trajectories were then compared to the actual drifter track for the 48-h segment.

The motion of a particle in a two-dimensional velocity field is described by:

$$ \frac{{{\text{d}}r}}{{{\text{d}}t}} = u\left( {t,r} \right) $$
(1)

Where r(x,y) is the position of the particle and u = (u,v) is the Eulerian velocity at position r and time t. The prediction of the particle trajectory was done by integrating Eq. (1) using a predictor-corrector scheme with the velocity given by:

$$ u = U + u\prime $$
(2)

Where U was the spatially interpolated radar velocity derived from either the UWLS or OI algorithms to the particle location and u′ is the contribution of the sub-grid scale variability and the measurement uncertainty. The model used to determine u′ in this study was based on a random flight model (Ullman et al. 2006; Griffa 1996). The random flight turbulent parameters that we used to determine u′ are shown in Table 1. The parameters are consistent with those calculated by Ullman et al. 2006 and used in the present application in SAROPS by the USCG. At each hour of each simulation, we calculated the location of the centroid of the particle cluster and the 95th percentile confidence region (Fig. 2). This confidence region was calculated by rank ordering the spatial bins, by frequency, of a two-dimensional frequency histogram of the final simulated particle locations (Ullman et al. 2006). In addition to the drifter trajectories, we also calculated the mean separation and 95th percentile separation of the actual drifter position relative to the initial drifter position for each 48-h segment. This position persistence scenario represents a search based on the last known position in the absence of surface current data. The processing and analysis metrics are consistent with those introduced in Ullman et al. (2006) and applied in Gong et al. (2010).

Table 1 Parameters used in the random flight model to determine particle trajectories
Full size table
Fig. 2
figure 2

Predicted particle locations (blue dots) using the OI currents 48 h after the start of the segment. The track of the centroid of the particles (green) and the actual drifter, SLDMB 43484 (red), is also shown. For reference, an instantaneous HF radar vector map is shown as overlaid black vectors

Full size image

4 Expected differences between HF radar and in situ data

HF radar provides a unique measure of the surface current. Before any evaluation can be done to assess the quality of HF radar surface current estimates, the expected differences due to environmental variability must be quantified (Graber et al. 1997; Kohut et al. 2006; Ohlmann et al. 2007; Paduan et al. 2006). The temporal averaging of the drifter and ADCP velocities was consistent with the HF radar processing to remove bias from different temporal sampling and averaging. Due to significant differences in the HF radar, drifter, and ADCP sampling in the horizontal and vertical, the same cannot be accomplished through spatial averaging.

4.1 Vertical differences

Since the HF radar currents are based on the signal scatter off the surface gravity wave field, the measurement is limited to the upper portion of the water column influenced by the waves. Depending on the assumed velocity shear profile, the effective depth of the 5-MHz system is approximately 2.4 m (Stewart and Joy 1974). The moored ADCPs sampled a bin approximately 2.0 to 3.2 ms deeper than this effective depth. The expected difference between the HF radar and ADCP observations will scale with the magnitude of the vertical shear between the effective depth and the depth of the ADCP measurement bin. Since the drifters are drogued to 1-m depth within the range of the HF radar measurement, the contribution of vertical shear to the observed differences is expected to be small.

The best measure of the vertical velocity shear that will contribute to observed differences between the ADCP surface bins and the HF radar observations at the time of this study were from the four moored ADCPs. For each ADCP, the velocity estimates from the bin closest to the surface were compared with the bin 3 m deeper. This depth range most closely matches the approximate distance from the surface HF radar measurement and the uppermost uncontaminated ADCP bin. This shallowest uncontaminated bin is defined as the shallowest bin without contamination from interaction with the surface. Given the ADCP frequencies and water depths, this bin was between 4.39 and 5.59 m below the surface for the four ADCPs. The root mean square (RMS) difference between these surface bins and the bins 3 m deeper was on average 2.1 cm/s in the east direction and 2.4 cm/s in the north direction. All values were between 1.9 and 2.9 cm/s with the square of the sample correlation coefficients (r 2) all greater than 0.96. The complex correlation indicates that not only are the velocities through the upper water column correlated (correlation coefficients of 0.97–0.99), but for each mooring the deeper current direction most correlated with the surface was offset by less than 1°. These comparisons indicate that the contribution of vertical shear to the measured difference between HF radar and the moored ADCPs is 2–3 cm/s.

4.2 Horizontal differences

Velocity shear in the horizontal contributes to the observed differences between HF radar and both the ADCP and drifter velocities. Again, the HF radar measurement is unique relative to the available in situ data in that it samples a spatial cell much larger than both the ADCP and individual drifters (Fig. 3). Given our range resolution of 5.85 km and azimuthal resolution of 5°, the spatial area of each radial velocity average scales with the range from the remote site. Ten kilometers from the site this area is approximately 7 km2 and grows to 54 km2 100 km from the site. In addition, both combination algorithms use the spatially averaged radial velocities within a radius on the order of 10 km in the total vector calculation. Shear in the flow across the scales sampled by the HF radar, drifters, and ADCPs will contribute to the observed differences between their current observations. These shears have been shown to vary significantly in time and space and contribute to the observed differences between HF radar and in situ sensors (Kohut et al. 2006; Ohlmann et al. 2007).

Fig. 3
figure 3

An expansion of the study site showing the location of the ADCPs (stars), the track of the February–March drifter (thick black lines), and the radial grid from the Loveladies, NJ HF radar coastal station (thin black lines). The triangle marks the deployment location of drifter #43484. For scale, the grid resolution is 5.8 km in range and 5° in azimuth

Full size image

Both the spatially separated array of ADCPs and clusters of drifters were used to quantify the horizontal shear in the surface velocity fields. There were three moorings deployed 10 km apart along a line perpendicular to the coast between the 40- and 60-m isobaths (Fig. 3). A fourth mooring was located 11 km upshelf along the same isobath as the center mooring in the cross-shelf line. The surfacemost bin of each mooring was compared to all other ADCPs to quantify horizontal differences across the scale of the array. For reference, the 10-km spacing is on the same order as the distance across adjacent range cells of the HF radar grid. The RMS differences between the surface bins of each ADCP were on the order of 5–7 cm/s in the east direction and 8–10 cm/s in the north direction. The complex correlation shows that there is a strong correlation (r 2 between 0.82 to 0.87) with angle offsets ranging from −5.9 to 4.4°.

The cluster of drifters deployed in April of 2007 was also used to provide an estimate of the horizontal shear. Unlike the ADCPs, the drifters move throughout the footprint of the HF radar network and the separation between them varies throughout the deployment. In the April deployment, the drifter separation ranged from 100 m to 100 km. At each 30-min time step, the current was estimated for each drifter in the cluster. From these time series of velocities, we calculated the RMS difference between drifters binned by range (Fig. 4). The RMS difference for a drifter separation of 10 km is based on all estimated velocity pairs from drifters between 0 and 10 km apart. As the distance between drifters increases so does the RMS difference. At the scale of the HF radar observation, the RMS difference in both the north and east velocity components is approximately 5 cm/s. This is consistent with the RMS difference determined from the ADCPs over a similar scale. As the distance between drifters grows toward 50 km, the RMS difference in each component reaches its maximum at 12–15 cm/s (Fig. 4).

Fig. 4
figure 4

RMS difference of the east (solid) and north (dashed) components of the estimated velocities for the drifters released in April 2007. All data between drifters within the distance plotted on the x-axis were used to determine the RMS difference

Full size image

The spatial contribution of known environmental variability is 5–10 cm/s based on the ADCP and drifter intercomparison. This gives us a measured quantity for the horizontal shears as they impact the evaluation and the application of each algorithm. These expected differences indicate that the environmental shears, particularly in the horizontal, can be a significant contributor to observed differences in the comparisons between HF radar and the available in situ currents and should be considered in this evaluation.

5 Results

5.1 UWLS and OI comparisons

The ADCPs were deployed in a region of greater than 90 % HF radar coverage with good remote site geometry. A scatter plot of the velocities between the center ADCP (MSF-C) compared to the UWLS and the OI base run 1 is shown in Fig. 5. The statistics between the surface bin of the four ADCPs and the UWLS vectors had RMS differences that ranged between 6.9 and 8.5 cm/s in the u component and 7.6 cm/s and 8.8 cm/s in the v component (Table 2). These statistics are meant to summarize the comparison between the in situ and remote currents and do not themselves reveal the episodic nature of the data mismatches. The data coverage for all comparisons is referenced to a complete in situ time series so that 100 % coverage means that all in situ observations have a concurrent HF radar observation. The data coverage for each ADCP/UWLS comparison time series was between 82.7 % for the offshore ADCP and 92.9 % for the inshore ADCP with at least 764 data pairs in each evaluation series. The OI run compared to the UWLS is run 1 with S x  = S y  = 10 and an error threshold of 0.95. For these OI total vectors, both the u and v components are consistent with the UWLS with slightly higher RMS differences from 7.6 to 8.8 cm/s and from 7.6 to 9.4 cm/s, respectively (Table 2, runs 1 and 5). All of these evaluations were based on at least 839 comparison pairs with a higher data coverage compared to the UWLS runs ranging between 98.3 % for the offshore ADCP and 99.6 % for the central ADCP.

Fig. 5
figure 5

Scatter plots of the UWLS (black) and OI (gray) compared to ADCP MSF-C for the east (top panel) and north components (bottom panel)

Full size image
Table 2 Summary of error statistics between HF radar and the moored current meters
Full size table

The drifter tracks crossed the 90 % coverage contours into regions of reduced HF radar coverage (Fig. 1). Scatter plots between drifter 43484 estimated currents and the UWLS and OI (run 1) currents are shown in Fig. 6. The comparison between this drifter and the UWLS totals was consistent with the ADCP evaluation with RMS differences of 7.2 and 9.7 cm/s in the u and v components, respectively. The coverage of the UWLS total vectors with the drifter, however, was significantly less with 418 data pairs (49.9 %). The OI total vectors compared to the same drifter had higher RMS differences of 9.3 and 12.6 cm/s in the u and v components, respectively. The coverage with the OI total vectors was significantly higher than the UWLS with 823 comparison pairs giving a data coverage of 98.3 %.

Fig. 6
figure 6

Scatter plots of the UWLS (black) and OI (gray) compared to drifter 43484 for the east (top panel) and north components (bottom panel)

Full size image

The April drifter tracks spent more time in regions of less than 90 % HF radar coverage. The spread in these tracks led to a significant variation in the comparison statistics. For the UWLS totals, the RMS difference ranged from 7.3 to 12.5 cm/s in the u component and 11.1 to 12.4 cm/s in the v component. The data coverage varied significantly between each drifter from 37.2 to 73.4 %. For the OI base run, the RMS differences were higher with a range of 9.4 to 12.4 cm/s in the u component and 11.6 to 14.4 cm/s in the v component. Like the February drifter, the percent coverage with the OI totals was significantly higher than the UWLS with a range across all April drifters between 55.7 and 96.7 %.

Drifter trajectories based on both the UWLS and OI base run solutions were compared to the actual track of drifter 43484. For each simulation up to 48 h, the distance between the centroid of the 1,000 simulated particles and the actual drifter location was calculated each hour (Fig. 7). The UWLS-based trajectories show a mean distance steadily increasing toward 12 km at 24 h and 18 km at 48 h after the start of each segment. The 95th percentile separation is a measure of the scale of the area encompassing 95 % of the released particles each hour. For the UWLS, this separation grows steadily with a maximum of 34 km at 48 h. The OI base run trajectories are shorter for all times with a mean distance of 9 km at 24 h and 14 km at 48 h. Similarly, the 95th percentile separation is consistently less than the UWLS solution with a maximum of 26 km at 48 h (Fig. 7). Both the UWLS and OI results fall below the drifter persistence values assuming a search based on the last known position alone.

Fig. 7
figure 7

Mean separation between the particle centroid and actual drifter (left) and the 95th percentile separation (right) for the UWLS (thick dashed) and OI (thick solid) simulations. For reference, we show the same statistics for the persistence simulation based on a static last known position (thin)

Full size image

5.2 OI sensitivities

Decorrelation scales

The decorrelation scale used to determine the weights for the input radials to the OI algorithm is defined as a distance east and north from the grid point. For our base run, this area was scaled as a circle with a radius of 10 km (S x  = S y  = 10). For the OI algorithm, a radial within and a little beyond this area are weighted based on an exponential decay with distance from the grid point. The rate of this decay is a function of S x and S y . As a variation to this base setting, we added asymmetry to the scale and oriented it to match the local topography. For this run, S y was rotated along the isobaths and stretched 2.5 times S x (S x  = 10, S y  = 25 and orientation = 31°CW from true north). For both the ADCP and drifter estimated currents, the comparison between runs 1 and 2 are very similar (Tables 2 and 3). The RMS difference and data coverage for all in situ observations are within 0.5 cm/s and 0.5 %, respectively (Tables 2 and 3, runs 1 and 2). Similarly, the drifter trajectories show very little sensitivity to the decorrelation scales. The mean separation between the centroid of the simulated particles and the 95th percentile separation follows similar trajectories to the base run with a separation of 9 km 24 h out and 14 km 48 h out (Fig. 8, top panels).

Table 3 Summary of error statistics between the HF radar and drifters
Full size table
Fig. 8
figure 8

Mean separation between the particle centroid and actual drifter (thick solid) and the 95th percentile separation (thin solid) for the OI simulations listed in Tables 2 and 3 for run 1 (upper left), run 2 (upper right), run 3 (lower left), and run 4 (lower right). For reference, the mean separation (thick dashed) and the 95th percentile separation (thin dashed) for the persistence simulation based on a static last known position are shown in each panel

Full size image

Error thresholds

For each of the OI runs described above, a lower threshold of 0.6 was applied in runs 3 and 4. We selected this lower threshold of 0.6 to optimize the quality of the data without significantly impacting coverage (Fig. 9). While there is an increase in the correlation with in situ observations, these lower thresholds come at a significant cost to data coverage. We selected 0.6 as a second error threshold to test the sensitivities because it has higher correlation than the 0.95 values and it falls just above the sharp decline in data coverage seen in values less than 0.5 (Fig. 9). The comparison with the ADCPs for these second set of runs showed a slight reduction in the RMS difference on the order of 0 to 0.5 cm/s with a nominal loss in coverage on the order of 5–8 % compared to runs 1 and 2 with a threshold of 0.95. As with the decorrelation scales, the error thresholds had a minimal impact on the drifter trajectory results (Fig. 8, bottom panels). In both runs 3 and 4, the lower threshold did not significantly change the mean of the distance between the actual drifter and the particle centroid or the 95th percentile separation seen in run 1. In each case, the mean gradually increases over the 48-h simulation to a centroid distance of about 14 km with similar 95th percentile separation.

Fig. 9
figure 9

Correlation (dashed) and data coverage (solid) for the comparison between the OI algorithm and drifter 43484 based on the normalized uncertainty of the total vector components. The threshold of 0.6 for runs 3 through 4 is shown as a gray line

Full size image

6 Discussion

6.1 UWLS and OI comparisons

The in situ ADCP array was centered in a region of greater than 90 % HF radar data coverage and good geometry. In this region, the comparison between the OI and UWLS algorithms was very similar in both RMS and percent coverage. All RMS difference and correlation (r 2) statistics were within 0.7 cm/s and 0.06, respectively with a slight increase in data coverage with the OI currents (Table 2, runs 1 and 5). Based on the ADCP comparison, the two algorithms can be interchanged with little impact on the resulting total vector quality. The location of the ADCPs in the center of the HF radar footprint ensured that there was consistent radial vector coverage at nearly orthogonal angles at every time step so that the total vector calculations will be based mostly on radials very close to the total vector grid point. Given these conditions, the mathematics associated with either algorithm gives very similar results. It is not until the comparison data move away from this optimal location that the strengths of each algorithm become more evident.

The drifter tracks took the in situ observations throughout the coverage, from regions above and below the 90 % contour and regions of good and bad radial site geometry. Across these regions of the HF radar footprint, the UWLS- and OI-derived velocity fields were more distinct. The UWLS algorithm had consistently lower RMS difference, typically 1–3 cm/s compared to the OI total vectors. While these OI-derived vectors were typically associated with lower uncertainty, the data coverage of the UWLS fields was significantly lower. The OI algorithm did fill the gaps seen around the edges of the coverage more effectively with data returns for the drifters on the order of 10–30 % greater than the UWLS total vectors.

The drifter estimated current comparisons highlight the distinction between the two algorithms. The UWLS, as it is configured for this study and consistent with the operation in the present national network implementation, will scale all vectors within the 10-km radius with an equal weight of 1. Any radial vector falling outside of that radius will have an equal weight of 0 and will therefore not be included in the current vector estimate. The OI implementation as tested here weights each radial vector based on its distance from the grid point. So a vector located 10 km away will have a weight of 0.37. All vectors inside of this distance will be weighted higher, and vectors that fall beyond will be weighted less. Unlike the UWLS implementation, these far field vectors will get included in the vector estimate, but will be weighted far less than the vectors that fall inside the decorrelation scale set by S x and S y . Given this, the UWLS implementation will only have a vector to compare to the drifter if there are radial velocities within 10 km of the grid point closest to the drifter. At those times, the OI will weight these velocities much higher than those that fall beyond 10 km. For the times when there is no UWLS solution but there is an OI solution, there cannot be radial velocities within 10 km of the closest grid point to the drifter. Therefore, the velocities that are included in the OI solution must fall outside the decorrelation scale of 10 km. The result is a total vector to compare with the drifter velocity that is based only on radial velocities more than 10 km from the grid point. The larger the horizontal shear between the grid point location and the input radial velocities, the greater the difference with the drifter velocity. If the shears are small over these scales, than a radial vector weighted relatively high, but further away from the grid point, will provide a representative velocity component for the estimated total vector. The 5–7-cm/s RMS difference observed across a similar scale in the ADCP array and April drifter cluster reinforces the significance of the observed shears in our study region and helps explain the issues that arise when weighting radials further away from the grid point too high.

As an example, we look at the comparison with drifter 43484. When we constrain the comparison to those times when there is a current estimate from the drifter, UWLS, and OI, we are limiting the HF radar solution to those times when there are radial velocities inside the 10-km radius set by the UWLS settings. In this case, the RMS difference between the UWLS and OI results compared to the drifter is 7.19 and 7.85 cm/s in the east direction and 9.70 and 9.12 cm/s in the north direction, respectively. In both the north and east directions, the OI results have come down to the lower uncertainties seen in the UWLS comparisons reported in Table 3. This shows a reduction in the RMS difference with the OI solutions of about 2–3 cm/s. So if there are radial vectors within the decorrelation scale of 10 km, the results from each algorithm are consistent. However, in those instances when there are no vectors within the 10-km scale, as is the case with the additional vectors seen in the higher data coverage of the original evaluation, the OI implementation will weight these far field vectors relatively higher and base the total vector calculation only on them. The result is a more filled in field with higher uncertainty.

The most striking difference between the two algorithms as implemented in this study was seen in the drifter trajectories. The UWLS simulations were consistent with those reported by Ullman et al. (2006) using data collected over the same area of the shelf 3 years before our study period. In both studies, the drifter mean separation was on the order of 10–12 km at the end of the 24-h segment. This is an improvement over the persistence solution in which we assume a search based on a static last known position (shown in Fig. 7). It is this improvement and the associated statistics that drove the inclusion of UWLS-derived total vectors into the first operational application of HF radar in SAROPS. It has been shown that this improvement with the UWLS can reduce search areas by a factor of 3 (Roarty et al. 2010).

For the OI simulations, the metrics for both the mean drifter separation with the particle centroid and the 95th percentile separation are better than those of the UWLS. The smaller 95th percentile separation across all hours of the simulation leads to a smaller search area, reducing search times and increasing the likelihood for a successful rescue. The OI-based trajectories have a mean separation of about 12 km 48 h into the simulation. This is of the same order as the UWLS results at 24 h into the simulation (Fig. 7). This improvement is based on the varying weights applied in the combination. The variable weight based on the distance from the grid point assures that radial velocities closer to the SLDMB location will be weighted more than those further away. The result is a more representative total vector to drive the virtual drifter trajectories. Based on both the intercomparison of the ADCPs and the drifter clusters, we quantify significant variation in the surface velocity over the decorrelation scales used in the OI. In the case of the UWLS, the vectors within 10 km of the actual drifter are all weighted the same, leading to radial velocities as far as 10 km away contributing equally to the resulting total vector calculation as those within 1 km. In the context of the observed shears, this will lead to a total vector less representative of the velocity at that grid point compared to the weighted OI results. The impact of weighted radials on the particle trajectories will scale based on the magnitude of the horizontal shears over the scales of the processing. Our objective in this study was to evaluate the two algorithms as they have been implemented in existing HF radar networks. There are of course opportunities to weight the radials vectors within a least squares approach. Based on this evaluation, this would likely move the weighted least squares estimated trajectories closer to the observed drifter track in a similar way as that seen in the OI results.

The UWLS solutions currently feeding the operational SAROPS tool for cases in the Mid-Atlantic Bight was shown to reduce search areas by a factor of 3 over other available data and forecasts (Roarty et al. 2010). These OI simulations indicate that 24 h from the initialization of these segments, the search area is reduced by an additional factor of 2. Smaller, more accurate particle trajectories led to smaller more accurate search areas that lower search times and increase likelihood of a successful search.

6.2 OI sensitivities

The UWLS algorithm has been applied to HF radar total vector calculations for at least 30 years. Since the implementation of the OI to calculate HF radar total vectors has only recently been introduced to the field, we tested some of the algorithm’s sensitivities to help determine the optimal application of the algorithm. The two categories of parameters we tested were the decorrelation scales (S x and S y ) and the threshold for the reported uncertainties in the east and north directions.

The impact of the selection of decorrelation scale was tested with both the ADCP and drifter estimated currents. S x  = S y  = 10 was our base run to be consistent with scales currently used in the UWLS implementation and a second parameter set (S x  = 10, S y  = 25) with an elongated scale oriented in the direction of the isobaths. This second parameter set incorporated the characteristics of the local dynamics into the calculation of the total vector. The ADCP comparisons showed very little sensitivity between the two runs (Table 2). In all the comparisons, the RMS differences were within 0.24 cm/s. Once again, the ADCPs are located in a region of HF radar data coverage greater than 90 %. In this region, the estimated totals will be based on vectors close to the grid point that are weighted much higher than those further away. So the sensitivity of the estimated total to the decorrelation scale is minimized. The fact that this high density region showed little sensitivity is a validation that the larger scale we chose was dynamically consistent. Care was taken in the selection of the parameters to define the asymmetry and orientation of the larger search area. Both the orientation and asymmetry were based on the dynamics of the region. The velocity fields of the MAB tend to have longer correlation scales along the isobaths compared to across the isobaths. The larger number of radial vectors included in the along isobath direction did not introduce more uncertainty into the estimated total vector because they still represent the velocity over the grid point.

The component-based normalized uncertainty reported in the OI algorithm is potentially very useful in not only assessing the quality of each total vector, but also in considering the component-based uncertainty in the application of the data. As part of the sensitivity tests, we evaluated the impact of different error thresholds on both the RMS difference and data coverage. While there was slight improvement seen in most of the runs, the sensitivity of the total vector field resulting from the higher or lower error threshold was minimal. In some cases, the RMS difference was slightly higher with the lower threshold. It is not until the thresholds fall below 0.5 that there is a sharper increase in correlation (Fig. 9). This of course comes at a much more significant reduction in data coverage.

The sensitivity of the different OI parameter sets was also small in the 48-h particle trajectory simulations. All four runs had a steady increase in distance between the actual drifter and the centroid of the simulated particles over time with a maximum of about 14 km at the end of the 48-h segments (Fig. 8). The larger decorrelation scale stretched in the direction of the isobaths slightly lengthened the region of the highest weighted radials around the drifter. This nominal change had little impact on the drifter trajectories. Consistent with the Eulerian results, the lower error threshold of 0.6 did not eliminate a large portion of the original data. Operators of regional and national networks must balance the improved RMS difference and reduced coverage to determine an appropriate error criterion (Fig. 9). At a value of 0.6, we saw some improved RMS difference with minimal impact on the data coverage. If a lower criterion were selected, we would expect a large drop-off in data coverage with a measured improvement in RMS difference.

7 Conclusions

Regional scale HF radar networks are developing into coordinated national systems to support a variety of coastal applications, from search and rescue and spill response to oceanographic research. For each application, it is important to understand and quantify the uncertainty of each total vector included in a given surface current map. Based on the evaluations conducted here, we identified the strengths of two existing algorithms to combine radial velocities to maps of total vector currents. Both the UWLS and OI algorithms produced surface current vectors that compared well with in situ sensors in the context of known sampling bias due to environmental variability. In regions of good geometry and greater than 90 % data coverage, they produced similar results with no significant difference between the two. The difference in the total vector solutions from the two algorithms was much more evident with the drifter comparisons where in situ observations moved throughout the HF radar data footprint. The largest contributing factor to the observed differences between the two algorithms is based on the known horizontal shears in the study region and how these particular algorithms weighted radial velocities. Given the significant shear, the evaluation highlighted the importance of weighting radials closer to the grid point higher than those that fall further away. It also showed the risk in estimating total vectors beyond the decorrelation scale in the absence of vectors closer to the grid point. In these cases, the OI algorithm filled gaps, but the interpolated data carried with it higher uncertainty. The most significant distinction between the two was seen in the drifter trajectories. With these simulations, the weighted radial approach of the OI was very effective at reducing both the distance between the actual drifter track and the scale of the search area. In this study, the OI further reduced the already improved UWLS-based search areas by an additional factor of 2. The results indicated that the OI output was relatively insensitive to the varying decorrelation scales and error thresholds tested here. In regions like the Mid-Atlantic Bight off the east coast of the USA with significant horizontal shears, we are able to quantify the importance of weighting radial velocities when combining into maps of total vectors. Operators of the large regional and national networks will need to decide on what processing algorithm and parameter sets should be used to ensure quality data output that will support activities like search and rescue and oil spill response.