This study investigates the results of blending altimetry-based surface currents in the Gulf of Mexico with available drifter observations. Here, subsets of trajectories obtained from the near-simultaneous deployment of about 300 Coastal Ocean Dynamics Experiment (CODE) surface drifters provide both input and control data. The fidelity of surface velocity fields are measured in the Lagrangian frame by a skill score that compares the separation between observed and hindcast trajectories to the observed absolute dispersion. Trajectories estimated from altimetry-based velocities provide satisfactory average results (skill score > 0.4) in large (~100 km) open-ocean structures. However, the distribution of skill score values within these structures is quite variable. In the DeSoto Canyon and on the shelf where smaller-scale structures are present, the overall altimeter skill score is typically reduced to less than 0.2. After 3 days, the dataset-averaged distance between hindcast and drifter trajectories, , is about 45 km—only slightly less than the average dispersion of the observations, km. Blending information from a subset of drifters via a variational method leads to significant improvements in all dynamical regimes. Skill scores typically increase to 0.8 with reduced to less than half of . Blending available drifter information with altimetry data restores velocity field variability at scales not directly sampled by the altimeter and introduces ageostrophic components that cannot be described by simple Ekman superposition. The proposed method provides a means to improve the fidelity of near-real-time synoptic estimates of ocean surface velocity fields by combining altimetric data with modest numbers of in situ drifter observations.
Accurate near-real-time estimates of ocean surface velocity fields are necessary for predicting upper-ocean biogeochemical transport and managing accident response efforts. This is especially true in the Gulf of Mexico (GoM), where highly developed fisheries and oceanic transportation routes coexist with intensive petroleum drilling efforts and tourism in a semienclosed sea subject to the frequent passage of tropical cyclones. Two massive oil spills, the explosion of the platforms IXTOC-I in 1979 (Jernelöv and Lindén 1981) and Deepwater Horizon (DWH) in 2010 (Crone and Tolstoy 2010), have occurred in Gulf waters. Frequent episodes of red tides and hypoxia have been induced by agricultural runoff of nutrient-enriched river water into the marine ecosystems (Sklar and Browder 1998).
At basin scales, the surface circulation in the GoM is mainly driven by the intrusion of the North Atlantic western boundary current from the Caribbean Sea (Morey et al. 2005). The warm anticyclonic inflow, called the Loop Current (LC), finds its way in the eastern GoM and displays a wide range of oscillations (Oey et al. 2005). Irregular shedding of Loop Current eddies (LCEs), their westward migration, and interaction with topography influence the mean anticyclonic flow (Cochrane 1972; DiMarco et al. 2005; Lipphardt et al. 2008). Eddy–shelf interaction is usually observed around the DeSoto Canyon, an erosional valley characterized by the right-angle intersection of the Mississippi–Alabama slope and the west Florida slope (Harbison 1968). In this region the eddy activity, the Mississippi River outflow (MRO), and occasional intrusions of LC/LCEs (Huh et al. 1981; Hamilton et al. 2000) induce an interplay between local and deep ocean flows, affecting cross-shelf transport (Vidal et al. 1992; Ohlmann et al. 2001; Morey et al. 2003; Hamilton and Lee 2005; Weisberg et al. 2005).
In operational situations, at moderate offshore distances, the primary data streams typically available for estimating upper-ocean velocity fields are altimetry data in the form of gridded composite fields from several satellites (Ducet et al. 2000; Le Traon and Dibarboure 2004; Bouffard et al. 2008; Dussurget et al. 2011; Rio et al. 2011; Vignudelli et al. 2011; Escudier et al. 2013); observed and modeled surface winds (Sienkiewicz and Ahn 2005; Chen et al. 2008; Plagge et al. 2008; Bricheno et al. 2013); and directed, drifter-based, in situ observations targeting local transport mechanisms (Price et al. 2003; Sharma et al. 2010; Breivik et al. 2013). Although complementary, satellite altimeter and drifter observations provide information about the surface velocity field at very different space and time scales. How to optimally combine these two data streams to produce composite surface velocity estimates for transport studies remains an open question.
Satellite altimetric data provides extensive spatial coverage and is capable of resolving large-scale to mesoscale structures with space and time scales of the order of 100 km and weeks. Presently, however, this global synoptic coverage comes at the expense of feature and dynamic resolution. Traditionally, surface velocity fields are obtained from altimetric data via geostrophy, implying that only the geostrophic component of the horizontal velocity field is captured (Wunsch and Stammer 1998). While the geostrophic balance holds for larger mesoscale features, ageostrophic contributions are increasingly significant at scales near and below the Rossby deformation radius (Capet et al. 2008; Klein et al. 2008). Because of the large spacing [O(100 km)] between satellite ground tracks (Ducet et al. 2000), submesoscale processes are not currently resolved by gridded satellite altimeter-derived sea level anomalies. Chavanne and Klein (2010) showed that even much higher-resolution (typically 6–7 km) along-track satellite data are subject to signal contamination from high-frequency motions, such as internal tides. The increasing interest and need for estimating surface advective transport at 10–100-km spatial scales over relatively short, days to weeks, time scales raises questions about the validity of using velocity estimates derived solely from satellite altimetry in this scale range. The Ekman component of the ageostrophic velocity calculated from wind stress forcing has been added to satellite altimetry velocity and tested in several global products (Lagerloef et al. 1999; Rio and Hernandez 2003; Sudre and Morrow 2008; Sudre et al. 2013). The resolution though is still limited by the wind forcing products (typically ¼° for satellite scatterometer observations and ~10 km for model outputs) and by the time scales of ocean response to winds (Sudre and Morrow 2008). Recent drifter-based observations in the DeSoto Canyon region clearly indicate the importance of local velocity fluctuations in setting dispersion rates at scales in the 100-m–100-km range (Poje et al. 2014).
In contrast to satellite-based altimetry, surface drifter observations provide direct estimates of the local surface velocity field. Coastal Ocean Dynamics Experiment (CODE) drifters are cross-shaped drogued buoys designed to follow sea surface currents within the first meter of depth (Davis 1985). With GPS tracking, finite position accuracy and errors in water-following capabilities produce velocity errors of 1–3 cm s−1 in moderate wind and wave fields (Poulain et al. 2009). Despite this accuracy, drifters only measure velocities along their trajectories. Drifter information is routinely used to infer statistical information on basin-scale velocity (Ohlmann et al. 2001; LaCasce and Ohlmann 2003; LaCasce 2008). In addition, drifter data have been used to improve altimetry-based estimates of geostrophic mesoscale velocities in boundary currents (Cuny et al. 2002; Centurioni et al. 2008), as well as to refine parameters for Ekman regression models used in global velocity products (Sudre and Morrow 2008). On synoptic scales, however, drifter-based reconstructions of surface velocity fields (influenced by both geostrophic and ageostrophic dynamics) has been hampered, even over modest spatial regions, by a lack of contemporaneous drifter measurements with adequate spatial data density.
In this paper we concentrate on hindcast estimates of the synoptic surface velocity field and particle trajectories in the eastern GoM during September 2012, approximately one month after the release of nearly 300 CODE drifters in the DeSoto Canyon region during the Grand Lagrangian Deployment (GLAD) experiment (Özgökmen 2012). Given the large number of drifters released over a short time period in a relatively small region of the ocean, the GLAD drifter dataset provides synoptic coverage of the surface ocean at various scales for nearly 6 months (Olascoaga et al. 2013). Direct comparisons between synthetic drifters advected by altimetry-derived velocities and the GLAD observations show visual agreement in overall mesoscale transport patterns from the canyon into deeper water (Olascoaga et al. 2013) but significant differences in Lagrangian dispersion statistics during the initial month after release (Poje et al. 2014). In the context of data-assimilating operational models, Jacobs et al. (2014) have used the drifter dataset to test basic assumptions in satellite data assimilation—in particular, background error variance amplitude and time correlations. By directly assimilating GLAD drifter velocities in a four-dimensional variational data assimilation approach (4DVAR), Carrier et al. (2014) and Muscarella et al. (2015) have quantified improvements in model velocity and trajectory estimates in the upper ocean.
Here GLAD drifter trajectories are blended with geostrophic velocities, as inferred by Olascoaga et al. (2013), using satellite altimetric sea surface height (SSH) data from Archiving, Validation, and Interpretation of Satellite Oceanographic Data (AVISO) subjected to a no-flow condition on the coastline. The results are assessed in terms of the fidelity of hindcast trajectories. The objective is to test a methodology that can be used in applications such as pollutant tracking or search and rescue activities where, in addition to available altimeter-based velocity fields, data from a limited number of directed drifter deployments are also available. In such operational situations, where accurate near-real-time trajectory estimates are required, an optimal blending of available drifter and altimeter observations provides direct data-only surface velocity field estimates while avoiding issues of model bias or systemic model error inherent in predictive estimation. Since such applications are focused on synoptic and regional scales, the data synthesis approach required is necessarily different from that used to combine altimeter and drifter observations to compute global products (Sudre and Morrow 2008; Maximenko et al. 2009).
Various methods have been proposed in the literature to reconstruct velocity fields from available trajectory information (Toner et al. 2001; Chang et al. 2011). Here we use the Lagrangian variational analysis (LAVA) approach (Taillandier et al. 2006a,b, 2008, 2010), which allows for statistically robust reconstructions of velocity fields either directly from purely Lagrangian observations or from combinations of Eulerian model/data and Lagrangian datasets. While LAVA has been previously applied to velocity fields from models and high-frequency (HF) radar (Taillandier et al. 2010; Chang et al. 2011; Berta et al. 2014), the proposed application presents a number of novel aspects. Blending drifters and satellite altimetry velocities is especially challenging because of the disparity in the spatiotemporal scales resolved by the two platforms. The extensive GLAD dataset allows for an unprecedented level of quantitative assessment, not only of the LAVA performance but also of the AVISO-based fields that are used as benchmark. Finally, a technical improvement of the LAVA method is presented that allows automated processing of spatially dense drifter data streams by clustering and averaging trajectory information when necessary.
When applying LAVA, the space and time scales used in the blending have to be chosen a priori. Here, we are interested in mesoscale variability that is expected to be potentially not well resolved by present altimeters. We focus on subinertial scales in time, filtering the data at 48 h, and we introduce a blending space scale R of the order of the estimated Rossby radius —that is, approximately 40 km in the open ocean and about 10 km in the DeSoto Canyon and shelf area (Chelton et al. 1998). The difference between the blending scales (10–40 km in space and longer than 1 h in time) and the satellite altimetry-resolved scales, of the order of 100 km in space and 1 week in time, suggests that the blending will allow for refined estimates of large mesoscale eddy variability in the open ocean, as well as significant modification of smaller structures in the DeSoto Canyon and shelf areas.
An important issue to be addressed is how to evaluate the results. In case of drifter assimilation (Fan et al. 2004; Lin et al. 2007), validation is often performed by first using assimilated drifters themselves and then considering other types of data, for instance, from subsurface acoustic Doppler current profiler (ADCP) measurements. In the case of blending, since the correction does not dynamically propagate and it is confined to the neighborhood of the observation, it is necessary to use data that are compatible with the blended ones and that are situated within the correction scale R. In our case, no other independent data of surface velocity (e.g., from HF radar or surface ADCP) in the area covered by the drifters were publicly available from the Gulf of Mexico Coastal Ocean Observing System data portal (http://data.gcoos.org). We therefore test the results with a subset of “control” drifters that are not used in the blending (Berta et al. 2014). The control drifters can be seen as pollutant proxies in operational applications, that is, substances carried by the currents whose position is not known and trajectories from the source need to be hindcasted. The main performance metrics are Lagrangian quantities in order to directly assess the quality of estimated hindcast trajectories. Additional Eulerian metrics are also used to characterize the changes induced in the satellite altimetric velocity field by the LAVA blending.
The paper is organized as follows: In section 2, the satellite, wind, and drifter data are presented. In section 3 the LAVA blending method and its GoM implementation are described; the trajectory hindcast calculations and the metrics used to evaluate them are defined. The results are presented in section 4, and conclusions and future perspectives are discussed in section 5.
a. Satellite data: AVISO-based fields
Several products such as AVISO (http://www.aviso.altimetry.fr; Rio and Hernandez 2004), Ocean Surface Current Analyses—Real Time (OSCAR; Johnson et al. 2007), and Geostrophic and Ekman Current Observatory (GEKCO; Sudre et al. 2013) are now available for global SSH and geostrophic velocities, based on multisatellite altimetric data. GoM surface velocities from these products have been tested by Liu et al. (2014) for trajectory hindcast using an 18-drifter dataset, and the results appear to be approximately the same for all products. Here we use the AVISO-based absolute geostrophic velocities as implemented in Olascoaga et al. (2013), with a spatial grid of about ⅜° and a time interval of 24 h. These fields are defined as the sum of (i) the mean dynamic topography (Rio and Hernandez 2004); (ii) the altimetric SSH anomaly distributed by AVISO; and (iii) a perturbation that guarantees that the normal projection of the velocity at the coastline vanishes (Iskandarani 2008), introduced in order to improve SSH in the nearshore region (Saraceno et al. 2008; Cipollini et al. 2010; Vignudelli et al. 2011).
In Fig. 1 the basic statistics from the satellite fields characterizing the circulation during the month of September 2012 are shown. The monthly mean of the following quantities are depicted: SSH anomaly (Fig. 1a), surface geostrophic velocity (Fig. 1b), SSH standard deviation (Fig. 1c), and SSH gradient magnitude (Fig. 1d). At the beginning of September, the northern boundary of the LC is found at ~24°N (Figs. 1a and 1b). This condition was already observed by Hamilton et al. (2005) and Schmitz (2005) after the LC extended northward, generally up to 26.5°–27°N, and an LCE detached from the LC (Sturges et al. 2005). Such an event occurred just before the GLAD experiment at the beginning of July 2012. After the shedding of an LCE, the penetration of the LC in the GoM may be further inhibited by the interaction with peripheral cyclones during the so-called blocking process (Zavala-Hidalgo et al. 2002). Figures 1a and 1b show the presence of cyclones just north of the LC, as well as the previously detached anticyclonic LCE in the central basin. We concentrate on the region covered by the drifters in the eastern GoM (Figs. 3a and 3b). A strong anticyclonic structure is evident (Figs. 1a and 1b) with a main LCE core around 25.5°N, 89°W and a smaller northwestern recirculation (~27°N, ~90°W). To the east of the LCE, a cyclonic region can be seen, with an intense southern eddy at 24°N, 86°W and an extended recirculation north of approximately 25°N reaching the MRO. This cyclonic structure is located just south of the 2500-m isobath, with the northern flowing branch approximately located at the southeastern margin of the DeSoto Canyon. The highest temporal variability (Fig. 1c) is found at the eastern and northern edges of the LCE (~22 cm at ~25°N, 87°W and ~27°N, 89°W), while the highest values of SSH gradient magnitude (Fig. 1d) correspond to the eastern and southern parts of the LCE and to the LC (~21°–25°N, 80°–88°W).
b. Wind data: The NCEP-NAM products and Ekman correction
Traditionally, the ocean’s surface velocity field has been approximated as a superposition of geostrophically derived and wind-driven components (Ekman 1905). Wind-driven Ekman currents result from the balance between the frictional stress due to the wind and the Coriolis force. The horizontal transport associated with Ekman currents has been found to significantly contribute to drifter trajectory patterns at 15-m depth (Lagerloef et al. 1999; Ralph and Niiler 1999; Lumpkin and Garzoli 2005).
where ϕ indicates the latitude, and and θ indicate the wind intensity and direction at 10-m height, respectively. Following Liu et al. (2014), this parameterization is applied to the altimeter data using wind fields supplied by the National Centers for Environmental Prediction (NCEP)’s North American Mesoscale Forecast System (NAM; NCEP 2014; Rogers et al. 2009).
NAM products have a spatial resolution of 12 km and a temporal resolution of 3 h. Surface Ekman currents are superimposed on the AVISO-based geostrophic velocities, and the new velocity field, denoted AVISO–NCEP, is used to evaluate the effect that the wind-driven component of currents has on Lagrangian transport estimates.
The average wind conditions during September 2012 (Fig. 2) are characterized by easterly winds (meteorological convention), which is the typical wind regime present during summer (Morey et al. 2005). This tropical weather pattern is occasionally influenced in summer by the rapid passage of weak cold fronts from the north. Higher wind variability is found in the northern part of the GoM as indicated by variance ellipses (Fig. 2).
c. Drifter data: The GLAD dataset
During the GLAD experiment (17–31 July 2012), approximately 300 CODE drifters were deployed and reported their GPS position every 5 min. The GLAD drifter dataset is publicly available online (Özgökmen 2012). CODE drifters are designed to closely follow currents within the first meter of the water column. Comparison with current meters shows that errors are within 1–3 cm s–1 for winds up to 10 m s−1 (Davis 1985; Poulain 1999; Poulain et al. 2009). To adequately sample the scales spanning the meso-/submesoscale transition, drifters were released according to a multiscale approach for which deployment sites were spaced at 2 km, with each site containing nine drifters arranged in triplets of nested equilateral triangles, with separations of 100 m between drifters within a triplet and of 500 m between triangles within a site. Deployment sites were chosen to cover the area of the DWH spill in the DeSoto Canyon. Further details of the GLAD deployment scheme, chosen to assess transport and dispersion in the range of 100 m–100 km, are found in Jacobs et al. (2014) and Poje et al. (2014). At the beginning of September, about 230 drifters were still reporting their position, with this number decreasing to nearly 170 by the end of the month. A map of the concentration of drifter data during September 2012 is shown in Fig. 3a, with a bin size of 0.25°.
The raw drifter data were treated to both remove outliers in position and velocity and also to fill occasional temporal gaps using a noncausal spline interpolation. The trajectories were low-pass filtered with a 1-h period cutoff and sampled at uniform 15-min intervals (Yaremchuk and Coelho 2014). For this specific application, the available dataset was further filtered to remove inertial oscillations (ranging from 24 h at ~30°N to 35 h at ~20°N; see Jarosz et al. 2007; Anderson and Sharma 2008) using a 48-h running mean. Trajectories were subsampled every hour in order to perform time integrations within the blending procedure. As further discussed below, for the LAVA application the complete drifter set is divided into two subsets (Fig. 3b): one group was used in the LAVA blending (b-drifters) and the remaining set was used as control data to quantify the effect of LAVA on transport estimates (c-drifters).
a. LAVA algorithm and implementation
LAVA is a variational algorithm used to blend Eulerian velocity fields with Lagrangian data represented by drifter trajectories. Here LAVA is applied to AVISO-based fields, described in section 2a, producing the blended fields indicated as AVISO-LAVA in the following.
The AVISO-based first-guess velocity fields are corrected by minimizing the distance (misfit) between observed drifter positions and numerical positions computed by advecting trajectories in the flow field. The correction is centered on the position of the drifter, and it is spread over a range R through finite iterations of the diffusion equation (Derber and Rosati 1989; Weaver and Courtier 2001). This procedure is implemented over successive time sequences .
The value of the parameters R and is dictated by the dynamics of the basin over which LAVA is applied and by the scale of the flow that is targeted. The space scale R usually corresponds to the Rossby radius in the area, while the time scale has to be shorter than the typical Lagrangian time scale of the drifters (Taillandier et al. 2006a). There are also two other operative parameters—that is, the grid size of the discretized velocity, which has to be smaller than R in order to resolve the features (), and the time step over which the data are provided, which has to be smaller than (). In Taillandier et al. (2006a), an extensive sensitivity analysis on the two main parameters R and has been performed, showing that results are robust for changes of R up to 50% and for .
The application of LAVA in the GoM requires the division of the whole area covered by drifters into two subregions characterized by different dynamics—southeast (SE) GoM and Mississippi, Alabama, and Florida (MAFLA)—as indicated in Fig. 3c. The difference in spatiotemporal scales between the two regions requires different choices of the LAVA parameters R and . The SE GoM, defined as a wide area in between 21°–27°N and 84°–92°W, covers the Yucatan Channel entrance and the Campeche Bank and is centered on an area of deep sea. On the other hand, the MAFLA area covers part of the shelf facing the Mississippi, Alabama, and Florida coastline, as well as part of the DeSoto Canyon (~29°N, ~87°W) and the Mississippi River delta (~29°N, ~89°W). The geographical limits span within 27°–30.5°N and 84°–91°W, so that the southern edge of the MAFLA region coincides with the northern border of the SE GoM area.
Chelton et al. (1998) estimate 40 km for GoM deep waters and 10 km for the shelf and slope area, and these values are used for the LAVA parameter R in the SE GoM and MAFLA. The grid size is chosen to be ⅜° in the SE GoM, corresponding to the grid size of the AVISO-based currents. In MAFLA, the AVISO-based velocity is linearly interpolated on a regular grid with resolution to allow adequate resolution of the smaller-scale shelf features. Given Lagrangian time scales 1–3 days (Ohlmann and Niiler 2005), the analysis time scale was set to 4 and 6 h for the MAFLA and SE GoM areas, respectively. In both cases the temporal resolution is given by the time step of low-pass-filtered drifter positions (1 h). The daily AVISO-based current maps are repeated hourly, as in operational applications the most recent velocity field is used until an updated map becomes available.
Because of the spatial inhomogeneity of the drifter data, the number of b-drifters available for the blending is different in the two selected regions (Fig. 3b). Moreover, the number of drifters in each region varies in time as drifters leave and enter the fixed domains. The average number of blended drifters () is 99 for the MAFLA area and 58 for the SE GoM, while the c-drifter subset is composed of 30 drifters in total. Control drifters represent approximately 15% of the GLAD drifters in September 2012 and are chosen to give an approximately homogeneous coverage of the eastern GoM. The average distance between c-drifters and b-drifters, , is the main parameter that characterizes the data coverage with respect to the target trajectories (Berta et al. 2014). A sensitivity study in Berta et al. (2014) showed that blending results indeed deteriorate at increasing but that errors are limited for . In this application, the average is approximately 14 km in the SE GoM area and 4 km in the MAFLA area, that is, smaller than , providing a test case that is expected to be effective and at the same time affordable in practical applications. The LAVA parameters for both applications are summarized in Table 1.
A visual example of the effects of the LAVA blending is shown in Fig. 4, where a comparison between the average AVISO field (Fig. 4a), the AVISO-LAVA blended one (Fig. 4b), and their difference (Fig. 4c) are shown. The spatial distribution of the effects of the blending depends on the drifter coverage during the selected days. The difference in the parameter R between MAFLA and SE GoM is evident in Fig. 4c, with the blending scale much more extended in SE GoM than in MAFLA. Differences between the AVISO and AVISO-LAVA blended fields, computed using the weighted average defined in section 3b, reach values of the same order of magnitude as the current itself, especially in the MAFLA region (Fig. 4a).
The LAVA algorithm is based on the assumption that the flow is characterized by a main scale of motion R that is resolved on a fixed grid. In reality, however, ocean flows are inherently multiscale. For our area of interest in the Gulf of Mexico, for instance, Poje et al. (2014) have shown that in addition to the mesoscale, there is significant submesoscale contribution to the overall dispersion. Because of flow variability, drifters within a given grid can have contrasting velocity information. This is a common problem in blending and assimilating data at high resolution and concentration, and in the case of Eulerian methodologies it is often treated simply by averaging the data in space and time (Dobricic et al. 2010; Poulain et al. 2012) or by grouping drifters according to their relative distance (Koszalka et al. 2011). Nevertheless, for a Lagrangian blending methodology, there is at present no standard approach. For small drifter datasets, trajectories can be manually selected, that is, chosen from far enough deployment sites so that drifters have a relative distance greater than about , in order to avoid conflicting velocity information at the grid scale in the blending process (Berta et al. 2014). On the other hand, for extensive datasets, such as GLAD, an automated procedure is necessary. Here we implement a simple method to perform averaging on clusters of trajectories, prescreening the drifters in order to maximize coverage while minimizing redundancy in the trajectories.
The procedure is based on two conceptual steps performed at each cluster average time . Here is chosen as = 2, to ensure that trajectory redundancy is entirely eliminated over the analysis period. The first step consists of identifying “clusters” (Lee et al. 2007; Pelekis et al. 2011), defined as an ensemble of trajectories that, during 2, maintain a separation smaller than some minimum distance, here defined as . All drifter positions belonging to a cluster are averaged into a single trajectory according to their center of mass.
The second step is motivated by the fact that for each 2 period, trajectories that do not belong to a cluster may still encounter, at discrete times , other drifters with separation less than . In these cases, for each encounter we select the trajectory with higher information content and discard the other. This is done by ranking the trajectories by information content, defined as the number of encounters , where corresponds to the maximum information possible. When trajectories with the same have an encounter, the selection is arbitrary. The total number of discarded trajectories is typically less than 5% of the whole GLAD dataset.
b. Performance metrics
We compute hindcast trajectories from the three different velocity fields: the AVISO-based fields (AVISO), the Ekman-corrected fields using NCEP-NAM (AVISO–NCEP), and the GLAD drifter blended fields (AVISO-LAVA). For each drifter trajectory, a numerical particle is initialized every 24 h at the observed position and integrated forward in time for 72 h. In all cases the trajectory computation is performed by integrating the Eulerian velocity field using a fourth-order Runge–Kutta scheme. The performance of each velocity field is evaluated using two metrics that compare numerical trajectories with in situ drifter trajectories.
Let us first indicate with D the separation between drifters and numerical trajectories, defined as
where are the components of the drifter position at time t, and the subscripts s and n indicate the in situ and numerical drifters, respectively.
We then indicate with the absolute dispersion of the drifters, defined as
where is the ratio of the separation between drifters and numerical trajectories and the absolute dispersion of the drifters after 72 h = 3 days. The 3-day period provides an (upper) estimate of the Lagrangian predictability time , and it has been chosen also in previous works (Ohlmann and Niiler 2005; Liu et al. 2014). The skill score is calculated for each drifter with numerical trajectories reinitialized at the observed drifter positions every 24 h. Along each observed trajectory, a skill score value is assigned every 24 h.
A second metric, , is given simply by the average, computed over all drifters at all times, of the separation between in situ and numerical trajectories in the 72-h period.
In addition to the Lagrangian metrics, we also compute Eulerian metrics to quantify the differences induced in the velocity fields by the LAVA blending. Even though the velocity fields are computed for each time step , in order to facilitate visual inspection of the results, averaged fields are considered by introducing the normalized average relative difference, , defined as
where and denote the AVISO and AVISO-LAVA surface velocities, respectively; and () indicates the average over the period p (area a). Two different average periods of days and days have been used in the MAFLA and SE GoM regions, respectively. This is due to the different typical Eulerian persistence time scales of the two regions: the MAFLA is influenced by weather synoptic variability of the order of a few days, especially in the slope and shelf areas (Weisberg et al. 2005), while the deep sea SE GoM is dominated by mesoscale eddies, which may persist up to a few months (Vukovich 2007). The velocity differences are normalized by the space- (a) and time- (p) averaged rms AVISO-based velocity in each region. The areas and periods over which the average is performed are limited by the drifter coverage (Berta et al. 2014). The same type of averaging procedure is also applied to the weighted average of the vectorial difference between AVISO and AVISO-LAVA fields (Figs. 4, 11, and 12).
In the following, the metrics described above (s and ) are presented for AVISO-, AVISO–NCEP-, and AVISO-LAVA-derived surface fields. In sections 4a and 4b, the complete GLAD dataset is used to benchmark the AVISO and AVISO–NCEP fields. In section 4c, where the AVISO-LAVA fields are considered, the GLAD dataset is partitioned into blended and control drifters.
Two complementary spatial maps of the skill score s metric for trajectory hindcasts obtained using the AVISO-based velocity fields are shown Fig. 5. The map of individual skill scores (Fig. 5a) demonstrates large spatial inhomogeneity of the drifters (Fig. 3a). As a consequence, it is difficult to accurately present individual skill scores over the entire region, especially in areas of high data density, where values are superimposed. In addition to the individual skill scores, s, binned average values, S (Fig. 5b), are computed using the same 0.25° bin size chosen for the drifter concentration map. To include also regions with low data concentration, no cutoff value or normalization on the number of data per bin is imposed (Fan et al. 2004; Liu et al. 2014).
The results in Figs. 5a and 5b are qualitatively similar and indicate the presence of clear gradients of skill score corresponding to different regions. The regions with the highest skill score (S up to 0.7–0.8) appear to be located in the strong eddies—that is, the LCE and the southern cyclone (see section 2a)—even though the coverage there is sparse. The strip between the two eddies characterized by a southward-flowing jet (approximately along 87°W and between 24° and 26°N) instead has low skills (). Another region with relatively high skill score () and with much higher coverage can be seen in the cyclonic region south of the DeSoto Canyon (~27°N, ~87°W). Conversely, regions with low skill score () are prevalent within the DeSoto Canyon and on the slope and shelf.
We expect that high skills correspond to the sampling of processes that are well resolved by satellite altimetry, that is, processes with a strong signal in terms of SSH and SSH gradient and with scales of the order of at least 100 km in space and 1 week in time. To investigate this hypothesis, in Fig. 6 we show separately the bins with high (; Fig. 6a) and low (; Fig. 6b) average skill score, superimposed to the monthly mean SSH gradient magnitude. A sensitivity study has been performed considering different skill score cutoff values in the range 0.3–0.7 and the results are qualitatively consistent. At first approximation, high skill score regions appear indeed to be correlated with persistent large mesoscale structures with high SSH gradient, such as the main eddies and the LC, while low skills areas are found mostly in smaller mesoscale and submesoscale regions, such as the interior of the DeSoto Canyon and the slope and shelf. At closer inspection, though, it appears that in some regions there is a significant variability, with a mixture of high and low skill bins. Examples are the southern cyclonic eddy (24°N, 86°W) and the strip between the anticyclonic and cyclonic regions as can also be seen directly from Fig. 5b.
The reasons for this variability are not completely understood at this time, but at least two mechanisms can be put forth. The first mechanism is related to the nature of dynamical processes. We can expect that within large mesoscale structures, and especially along their fronts, instabilities can occur with significantly shorter space and time scales with respect to the eddies themselves (Zhong and Bracco 2013). These processes can be characterized by ageostrophic velocities, and they are not correctly captured by satellite altimetry, so that the associated skill scores are low. The second mechanism is related to the characteristics of the observing system. Satellite altimetric coverage varies in space and time, and we can expect that periods of low coverage in our region of interest would correspond to lower skill scores.
We conclude the analysis by considering the average metric. The results are shown in Fig. 7 (red line), together with the (black dashed) line for comparison. The metric is slightly smaller than , but the difference is certainly not significant, given the size of the variability. This result suggests that, even though the skill metric is relatively high in certain regions, the overall distance between the hindcast and observed trajectories is very close to the average distance traveled by the drifters. This means that on average the improvement of using satellite-altimeter-derived trajectories is marginal with respect to using the zero a priori knowledge that assumes that particles do not move from their initial positions. Technically, the difference between the results in terms of and s is mainly due to s being set to zero anytime the distance between the hindcast and observed trajectories is greater than the travel length (i.e., no negative skill values are considered). Conceptually, the two metrics highlight different aspects. The skill score s allows for identifying the regions where indeed there is an advantage in using the hindcast trajectories, but it does not quantify the error that is made when the skill is null. The metric, on the other hand, provides bulk information on the average performance of the hindcast, while it does not provide information on regional differences. Each metric has its advantage and disadvantage and it is useful to characterize the results with both of them.
To evaluate the effect that the wind-driven component of the currents has on Lagrangian transport estimates, surface Ekman currents (estimated from the NCEP-NAM wind model) are superimposed on the AVISO-based geostrophic velocities. The map of binned skill score S in the AVISO–NCEP case (Fig. 8a) is qualitatively similar to the AVISO case (Fig. 5b) even though it shows some improvements in certain bins, especially in the shelf area. The slight enhancement is in agreement with the results by Liu et al. (2014). Nevertheless, a close look at the skill differences between AVISO–NCEP and AVISO (Fig. 8b) shows that in some cases the addition of the Ekman effect can also lead to lower skill score values, even though the net value is slightly positive. Similarly, the metric for AVISO–NCEP (Fig. 7, blue line) shows only very marginal improvement with respect to the AVISO case. In fact, the average separation between synthetic and real particles is about 45 km after 72 h, approximately the same distance as for the AVISO case and for the average absolute dispersion. The standard deviation for the AVISO–NCEP case spans the same range as for the AVISO case. Therefore, for this application, the addition of the Ekman effect does not significantly decrease the uncertainty of the Lagrangian transport.
Different possible concurrent reasons can be given to explain this result. First of all, the information contained in the AVISO-based currents and NCEP winds resolves scales on the order of 100 km and 10 km, respectively. On the contrary, drifters are likely to be influenced also by very localized forcings. Also, the open-sea area presents dominant geostrophic dynamics (LC and its eddies) (Sudre and Morrow 2008) and an Ekman component addition is not expected to be significant in the absence of strong frontal passages or hurricanes. It should be noted that winds were moderate during the examined period (Fig. 2). On the other hand, on the shelf and DeSoto Canyon, where the action of the wind is potentially more significant, the superposition of the Ekman component on geostrophic currents does not take into account the complex response to changes in wind forcing in terms of time scales (Stewart 2008; Sudre and Morrow 2008) as well as several other processes contributing to the surface current dynamics, such as the eddy-induced shelf break and slope circulation (Ohlmann et al. 2001; Wang et al. 2003; Hamilton and Lee 2005), river discharges (mainly Mississippi and Apalachicola) (Schiller et al. 2011; Kourafalou and Androulidakis 2013), upwelling events (Nowling et al. 2000; Hsueh and Golubev 2002), wind-driven currents from the west Florida shelf (Yuan 2002; Clarke and VanGorder 2013), and the submesoscale-induced transport (Poje et al. 2014).
In Fig. 9a, the concentration of the c-drifters over both the MAFLA and SE GoM regions is shown. As for the complete GLAD dataset (Fig. 3a), the highest concentration of positions is found close to the deployment area. The c-drifters cover most of the GLAD region, except for the LCE, where the original coverage was already sparse. Binned skill score values for AVISO and AVISO–NCEP from the c-drifters are shown for comparison in Figs. 9b and 9c, and they appear qualitatively similar to the complete results in Figs. 5b and 8a. Results from AVISO-LAVA are shown in Fig. 9d, and it is immediately evident that the LAVA blending significantly improves the performance. High skill scores are noticeable in both the SE GoM and MAFLA regions, including the cyclonic structure in front of the DeSoto Canyon, the strip within the southward jet between the anticyclonic and cyclonic eddies, as well as the DeSoto and slope and shelf area. The only area that is only marginally improved is the southern cyclone, characterized by high skills also in the AVISO case. Skill scores values for AVISO-LAVA are frequently higher than 0.8, and only a few bins have values lower than 0.4.
These results confirm previous outcomes obtained by applying the LAVA blending to velocity fields from models and HF radars (Chang et al. 2011; Berta et al. 2014). Drifters directly sample transport by currents at various scales within the first meter of water depth, which is influenced by very complex dynamics induced by air–sea interactions, dynamical instabilities, and interactions with the MRO. Drifter blending, therefore, has the potential of complementing satellite altimetry fields at scales that are not sufficiently resolved, while refining resolved structures by introducing information on environmental variability as well as possible ageostrophic components.
It is interesting to look separately at the plot for the two areas, MAFLA and SE GoM, over which LAVA has been applied (Fig. 10). The two areas are characterized by very different spatiotemporal dynamical scales (Chelton et al. 1998; Leben 2005; Weisberg et al. 2005), and therefore we expect different trends for both and . In fact, the average distance traveled by drifters is about 60 km in SE GoM, while it is almost halved (approximately 36 km) in MAFLA. This difference is due to the fact that typical velocities in the shelf and slope region are lower than in the open ocean (Oey et al. 2005; Ohlmann and Niiler 2005). Velocities in the MAFLA (SE GoM) region are on average about 0.15 m s−1 (~0.3 m s−1), reaching 1 m s–1 within LC and mesoscale cyclones. The curve for the AVISO and AVISO–NCEP cases lies close to the line of (~52–55 km for SE GoM and ~40 km for MAFLA). Note that in this case, considering the reduced c-drifters dataset, the AVISO–NCEP is actually slightly higher than AVISO, even though the difference cannot be considered significant given the variability. The AVISO-LAVA curve shows significant improvements with a final average separation much lower than average absolute dispersion (~21 km for SE GoM and ~17 km for MAFLA).
The effects of LAVA blending on the AVISO velocity fields are illustrated for the SE GoM and MAFLA in Figs. 11 and 12, respectively. The visualization and the metrics are different from Fig. 4 because the two regions are shown separately to provide more details and also are averaged over different time periods, reflecting the typical persistency of dynamical structures in each area (Weisberg et al. 2005; Vukovich 2007). For the SE GoM a longer time averaging is used (15 days) with respect to MAFLA (3 days).
The SE GoM results (Fig. 11) show the average circulation in the second half of the month (16–30 September), when many drifters (Fig. 11b) moved southward following the jet between the cyclonic and anticyclonic eddies, and some of them got trapped in the southern cyclone, whereas other ones drifted northwestward following the anticyclone. The cyclone–anticyclone system is reproduced by the AVISO velocity field (Fig. 11a), but the LAVA blending induces significant differences, especially in the jet area. The weighted average intensity of AVISO currents, normalization term in definition [Eq. (5)], is about 0.35 m s−1. In Fig. 11c, reaches almost 200% in the area of the southward jet, whereas along the western margin of the LCE the difference can be locally of the order of 100%. For the remaining covered areas is lower, mostly below 60%, corresponding to a magnitude of ~0.21–0.28 m s−1. The vectorial velocity difference in Fig. 11d shows that LAVA blending significantly modulates the jet, inducing a more extended longitudinal shear, and impacts the two eddies even though at a lesser extent. In summary, AVISO appears to capture the large mesoscale structures but their details are introduced by the drifters. This is the reason for the great change in skill score between AVISO and AVISO-LAVA (Fig. 9), especially in the southern jet (from less than 0.4 for AVISO to 0.7–0.8 for AVISO-LAVA).
The MAFLA circulation during 22–24 September (Fig. 12) shows the presence of the northwestward flow south of the DeSoto Canyon in both the AVISO velocity (Fig. 12a) and the drifter trajectories (Fig. 12b). The circulation in the canyon and on the slope and shelf generally appears to be anticyclonic and quite complex, with marked differences between AVISO and the drifters. The drifters also suggest the presence of some smaller-scale features, such as local recirculations on the two sides of the Mississippi River (MR) delta (~29°N, ~89°W), and on the eastern DeSoto Canyon slope around 29°N, 87°W. The AVISO-LAVA field (Fig. 12c) shows significant differences with respect to AVISO, especially regarding the anticyclonic area. The weighted average intensity of AVISO currents, normalization term in definition [Eq. (5)], is about 0.19 m s−1. Values of reach 200% along the eastern side of the anticyclonic pattern and around the MR delta close to the shelf edge, where differences are on the order of 0.4 m s−1. Only in the northwestward flow are the differences are relatively small, less than 60%. In several areas, the vectorial velocity difference (Fig. 12d) is in the opposite direction and of the same order of magnitude with respect to the AVISO velocity, especially along the shelf break and along ~85°–86°W. This indicates that the AVISO field does not reproduce the smaller mesoscale structures of the DeSoto Canyon and of the slope and shelf, and that drifter blending induces extended changes in the velocity patterns. This is in agreement with the skill score results in Fig. 9. The northwestward flow south of the DeSoto Canyon is well resolved by the satellite altimeter and accordingly it displays small values of and high values of skill score. The canyon and shelf areas have low skill scores for AVISO (less than 0.4), whereas S increases to values generally higher than 0.6 up to 0.8 for AVISO-LAVA. Only the very few S bins on the northern shelf (Fig. 9) still have low skill score values, because the blended drifters coverage is very low.
5. Summary and discussion
The performance of trajectory hindcasts is evaluated against drifter trajectories observed during the GLAD experiment. We consider three velocity fields. The first two fields, similar to those considered by Liu et al. (2014), are AVISO-based geostrophic velocities and the same fields with the addition of an Ekman component from the NCEP-NAM winds, named AVISO and AVISO–NCEP, respectively. The third velocity field (AVISO-LAVA) is computed by the variational blending of AVISO data with a subset of GLAD drifter observations using the LAVA technique.
The first novel aspect here is the application of LAVA to satellite-altimetry-derived velocity fields. The second is the ability to blend large-scale altimetric fields with readily available, but highly localized, drifter data. Approximately one month after deployment, the GLAD trajectory dataset provides information from the submesoscale-rich DeSoto Canyon to the mesoscale-driven open ocean. As such, the performance of the data blending approach can be estimated across very different dynamical regimes. The large number of observations permits partitioning of the data into subsets for both input to the LAVA blending and control observations for performance evaluation.
The results are analyzed using two Lagrangian metrics: the nondimensional skill score s based on the normalized separation between individual hindcast and drifter trajectories over three days, and the time-dependent average distance, , computed over all the drifters in a given region. Eulerian metrics are also computed to evaluate the differences between the AVISO and AVISO-LAVA velocities due to the blending of trajectory observations.
Results for the AVISO-based fields show that the binned average skill score S tends to be higher () in open-ocean large structures that are well resolved by the altimeter, that is, characterized by high SSH and SSH gradient magnitude and with space and time scales of the order of 100 km and a week, respectively. This is consistent with the analysis based on Lagrangian coherent structures from AVISO-based velocity by Olascoaga et al. (2013). Regions characterized by less energetic and smaller mesoscale and submesoscale features, such as the DeSoto Canyon and the shelf, have typically reduced skill score using AVISO-based fields. This is in agreement with previous results by Liu et al. (2014). The high coverage provided by GLAD drifters, though, also shows that the variability in skill score is very high even in the open ocean and that high SSH gradients can correspond to low skill scores. In particular, the jet between the two main cyclonic and anticyclonic eddies is characterized by low skill scores, less than 0.4. This variability can be due to a number of reasons. On the one hand, dynamical processes can lead to the occurrence of velocity variability within the mesoscale structures that is not resolved by satellite altimetry. Examples are high horizontal shears, or instabilities with smaller space and time scales. On the other hand, more structural reasons related to the observational platform can also play a role. Satellite altimetry coverage varies significantly in time, and this can influence the results. The metric computed over the whole dataset shows that the distance between the hindcast and drifter trajectories is on average approximately 45 km, slightly smaller than the average distance traveled by the drifters, .
Results from AVISO–NCEP are similar to AVISO in terms of skill score and . Statistics on the complete dataset shows a small improvement over shelf areas, as in Liu et al. (2014), but it is not significant given the high variability. The physical reason for this result is most likely due to the fact that in our region of interest the dynamics are mostly influenced by mesoscale and/or submesoscale processes (Poje et al. 2014), for which wind action cannot be simply described as a superposition between geostrophic and Ekman flow (Nowling et al. 2000; Hsueh and Golubev 2002; Hamilton and Lee 2005; Clarke and VanGorder 2013; Kourafalou and Androulidakis 2013).
Finally, the AVISO-LAVA results show a significant improvement of the skill score in all dynamical regions, that is, in the open ocean as well as in the DeSoto Canyon and slope and shelf area. Skill scores are frequently higher than 0.8, and only a few have values less than 0.4. The values are of the order of 20 km with an uncertainty decrease of about 50% with respect to . An analysis of the velocity fields from AVISO-LAVA shows significant changes with respect to the AVISO velocity. Local differences between AVISO and AVISO-LAVA can approach 200% of typical velocities in both the open ocean and the DeSoto Canyon and shelf regions. The nature of the difference, though, varies according to the dynamical region considered. In the open ocean, the large mesoscale field estimated by AVISO is qualitatively consistent with AVISO-LAVA, but the blending introduces important modifications on the velocity structures. In particular the jet between the two main cyclonic and anticyclonic eddies is highly impacted by LAVA blending that introduces a more extended longitudinal shear. In the DeSoto Canyon and slope area, LAVA blending substantially modifies the velocity field, even changing velocity direction in some points, and introducing smaller structures that are not present in AVISO. This is consistent with the fact that in shelf areas dynamical scales are smaller and not adequately sampled by AVISO as in deeper waters. Drifter information, therefore, allows for reintroducing the high environmental variability of the near-surface (upper 1 m) circulation, including also the complex forcing interaction that is not described by the classical Ekman response to large-scale winds. This local variability may be undersampled by the satellite altimeter which, in return, provides large-scale features of the deeper circulation.
Looking at the dispersion plots from a different angle, useful considerations can be inferred concerning applications in the scenario of an accident at sea. Let us consider the (control) c-drifters as a proxy for the advected pollutant so that their absolute dispersion represents the distance contaminant particles have traveled from the source over a certain period of time. Thus, the average absolute dispersion measures the maximum uncertainty on particle positions. Consider now the case when AVISO (or AVISO–NCEP) currents are known and used to nowcast the pollutant patch by advecting synthetic particles from the contaminant source. The average separation between numerical trajectories and c-drifters (pollutant proxies) compared with tells us that the velocity information from AVISO (or AVISO–NCEP) acts to reduce the search range by approximately 8%–13% in the SE GoM. For the MAFLA area, AVISO (or AVISO–NCEP) currents do not improve the Lagrangian transport estimates. On the other hand, if we consider the Lagrangian transport using LAVA blended fields, the values of for AVISO-LAVA suggest that the uncertainty in pollutant position decreases drastically, with the contaminant search range reduced by approximately 65% and 53% in the SE GoM and MAFLA regions, respectively. We also recall that the average distance between blended and control drifters, , is approximately 14 km in SE GoM and 4 km in MAFLA area, which is about half of the parameter for both experiments. This has important consequences when dealing with a real emergency scenario in which the exact position of the pollutant source is not known and mitigation procedures take place some hours after the accident so that drifters are typically launched some kilometers away from the actual contaminant position. Even in these cases, LAVA blending still provides considerable improvements of the Lagrangian transport estimates in the accident area.
In summary, the results confirm that trajectory hindcasts in the GoM open-ocean energetic mesoscale regions can be in first approximation satisfactorily estimated by satellite-derived fields. This is remarkable, since it indicates that large-scale geostrophic velocities can control the flow in the upper meter, which is subject to many complex processes. On the other hand, even within the mesoscale, the space and time variability cannot be resolved by satellites, and regions with smaller scales like the DeSoto Canyon and shelf have very limited altimetric skill scores. Drifter blending is a very effective way to complement satellite altimetric fields. The present results indicate that an affordable launching resolution of the order of half Rossby radius in the area of interest can be effective (see also Berta et al. 2014). The LAVA blending method has been demonstrated to be easily adaptable to any region, provided that the dominant dynamical scales are known, and therefore it is expected to be faster and simpler to implement than a full assimilation procedure.
This research paper, developed in the framework of CARTHE, was made possible by the grant from BP/Gulf of Mexico Research Initiative. This work was also supported by the international collaboration with the CNR-ISMAR, through the EU-MED project Tracking Oil Spills and Coastal Awareness Network (TOSCA). The authors are in debt to E. Coelho, M. Iskandarani, G. Novelli, A. J. Mariano, E. H. Ryan, F. J. Beron-Vera, G. A. Jacobs, B. K. Haus, H. S. Huntley, B. L. Lipphardt Jr., and A. D. Kirwan Jr., all of whom made the GLAD experiment possible. The authors wish to express their gratitude to Dr. Vincent Taillandier, developer and user of the first LAVA version, from which this new LAVA application has originated. The authors are grateful to Prof. Vassiliki Kourafalou and Dr. Matthieu Le Hénaff for the valuable discussion about GoM dynamics and satellite altimetry. Acknowledgments are due to Dr. Simone Marini for sharing his experience on clustering techniques. The authors thank Dr. Jean Mensa for communicating useful information about the SSH dataset.