Submesoscale fronts arising from mesoscale stirring are ubiquitous in the ocean and have a strong impact on upper-ocean dynamics. This work presents a method for optimizing the sampling of ocean fronts with autonomous vehicles at meso- and submesoscales, based on a combination of numerical forecast and autonomous planning. This method uses a 48-h forecast from a real-time high-resolution data-assimilative primitive equation ocean model, feature detection techniques, and a planner that controls the observing platform. The method is tested in Monterey Bay, off the coast of California, during a 9-day experiment focused on sampling subsurface thermohaline-compensated structures using a Seaglider as the ocean observing platform. Based on model estimations, the sampling “gain,” defined as the magnitude of isopycnal tracer variability sampled, is 50% larger in the feature-chasing case with respect to a non-feature-tracking scenario. The ability of the model to reproduce, in space and time, thermohaline submesoscale features is evaluated by quantitatively comparing the model and glider results. The model reproduces the vertical (~50–200 m thick) and lateral (~5–20 km) scales of subsurface subducting fronts and near-bottom features observed in the glider data. The differences between model and glider data are, in part, attributed to the selected glider optimal interpolation parameters and to uncertainties in the forecasting of the location of the structures. This method can be exported to any place in the ocean where high-resolution data-assimilative model output is available, and it allows for the incorporation of multiple observing platforms.
Submesoscale structures, with horizontal scales of ~1–10 km and vertical scales of ~100 m, are ubiquitous in the ocean (Thomas et al. 2008; McWilliams 2016). However, their spatial distribution is not homogeneous. Strong, highly dense submesoscales are found in regions where the strain field exceeds the relative vorticity field, for example, at the edges of mesoscale eddies and in elevated strain areas between eddies, but not inside the eddies, where the strength of the relative vorticity exceeds that of the strain field. Here we present a new strategy for identifying and autonomously sampling those regions where there is strong production of submesoscales (Fig. 1). As we will see in the next paragraphs, this strategy is different from other methods that focus on following (tracking) one specific front.
In the framework of the Satellites to Seafloor project (Thompson et al. 2017) funded by the Keck Institute for Space Studies (KISS; http://www.kiss.caltech.edu/new_website/techdev/seafloor/seafloor.html), we use data-assimilative primitive equation ocean model forecasts of currents to guide an underwater glider remotely controlled by a shore-based path planner. The glider path is autonomously chosen by decision-making software that extracts information on ocean features from the numerical ocean model, evaluates different potential glider paths in real time, and selects the most promising path based on “feature scoring” (defined in section 3b) (Fig. 1). This unique closed-loop science framework driven by an ocean model represents an advance in autonomous ocean exploration and monitoring. By analyzing all the possible virtual glider paths for a given day, and selecting the optimal choice, our method takes advantage of the model forecast to detect which part of the ocean is more relevant to sample (according to our defined features). This approach could have significance for the exploration of poorly understood regions of the ocean, such as under ice shelf or deep ocean surveys, as well as for coordinating a large global array—for example, hundreds—of piloted autonomous vehicles.
Autonomous vehicles, and in particular ocean gliders (Eriksen et al. 2001; Sherman et al. 2001; Webb et al. 2001), with their longer endurance and rapidly improving reliability (Davis et al. 2002; Brito et al. 2014; Rudnick et al. 2016), have been increasingly used to study the ocean. However, to optimize sampling we need to use tools that include a comprehensive view of the oceanographic setting (common in ocean observing systems) and real-time decision-making. Because of the relatively fast time scales at which the ocean submesoscale evolves, a method for optimally navigating gliders or autonomous underwater vehicles (AUVs) should ideally be fully automated, with minimum or zero human intervention.
Autonomous sampling of the ocean has experienced major advances since the Autonomous Ocean Sampling Network field experiment (Curtin et al. 1993; Curtin and Bellingham 2009; Ramp et al. 2009). The Adaptive Sampling and Prediction field experiment in Monterey Bay, off the coast of California (Leonard et al. 2010), is an example of coordinated control of multiple autonomous vehicles that consisted of adaptive, coordinated, model-based sampling with human decision-making. Coordinated control of multiple underwater vehicles has the objective of avoiding redundant measurements (Zhang et al. 2007). Adaptive sampling requires the use of observations or data-assimilative models to plan optimal trajectories (Curtin and Bellingham 2009). A necessary component for full autonomous ocean sampling is the use of a path planner to control the underwater assets (Fiorelli et al. 2004).
Several prior studies on autonomy focus on onboard decision-making to modify the vehicle path. For example, Cruz and Matos (2010) control the vertical position of an AUV based on recognition of entering and exiting the thermocline. Zhang et al. (2011) used onboard autonomy on an AUV to sample an oil spill. On one cast (up or down), the AUV detected a vertical oil feature (based on a peak hydrocarbon signal) and on the next cast it triggered a water sampler at that depth. In Zhang et al. (2012) the AUV analyzed the time evolution of the vertical temperature gradient to detect an upwelling front (a surface thermal front with isotherms slanted in the vertical). The boundary between stratified and homogeneous temperature regions defined the edge of the upwelling. The AUV then used the horizontal location of the boundary of the upwelling to define a lateral (x, y) zigzagging behavior to better sample the front. A similar technique was also used to determine the vertical location of the thermocline (Zhang et al. 2010), enabling subsequent vertical AUV casts to focus on this depth range. Cruz and Matos (2014) estimate the curvature of a horizontal front and then determine the lateral trajectory to cross it perpendicularly to optimize measurement of the cross-front structure. Sun et al. (2016) use an adaptive gradient method for thermocline detection and compare against peak gradient, average gradient, and human labeling. These techniques can be combined to focus on fronts in multiple dimensions (Zhang et al. 2016) and target both physical and biological processes (Godin et al. 2011).
Other work has used machine learning methods, in the form of policy learning (Magazzeni et al. 2014), to track features. This approach, however, requires a large amount of data or simulation. Other marine autonomy work focuses on flying formations relative to a tracked feature as indicated by a drifter (Das et al. 2012). Petillo et al. (2012) employs a hierarchical autonomy approach in which a model simulation of a plume is used to estimate the location of the plume boundaries and then vehicles are sent, with individual paths, to track these boundaries.
The objective of the Monterey Bay test experiment is to show that a data-assimilative forecasting numerical model can be used, in an autonomous fashion, to optimize and control an underwater vehicle sampling strategy. The outline of the paper is as follows. Section 2 briefly describes the oceanographic setting within Monterey Bay, defines the target features (thermohaline fronts), and introduces the Seaglider used in the ocean deployment. Section 3 provides a description of the methodology, which includes the Monterey Bay Coastal Ocean modeling system, the feature extraction method, the autonomous planner used to control the Seaglider path, and the fundamentals of the feature detection planner and path reformulation. The evaluation of the method is presented in section 4. The comparison of the glider and model results are discussed in section 5. In section 6 we discuss our work in the context of previous efforts performed in Monterey Bay. Final remarks are summarized in section 7.
2. Oceanographic setting
a. Monterey Bay dynamics
The Monterey Bay circulation is governed by the California Current System (CCS) (Hickey 1979; Lynn and Simpson 1987; Checkley and Barth 2009). The CCS is characterized by persistent coastal upwelling in response to prevalent northerly winds, which generate highly productive cold coastal areas. Offshore (>150 km), the California Current (CC) flows southward with surface speeds of ~0.25 m s−1. Near the coast (<50 km), the surface flow varies seasonally, flowing northward in fall and winter (Reid and Schwartzlose 1962), as the so-called inshore countercurrent (IC) (Lynn and Simpson 1987). The IC is intermittent in space and time, and it often receives regional names, such as the Davidson Current (Reid and Schwartzlose 1962) north of Point Conception (34.45°N, 120.47°W).
Below the IC, the subsurface California Undercurrent (CU) flows northward. South of Monterey Bay, at Point Sur (36.31°N, 121.90°W), the CU separates from the coast as a result of topographic curvature and flow inertia (Molemaker et al. 2015) and forms long-lasting mesoscale anticyclonic eddies (Fig. 2). Along the edge of these eddies, and likely caused by the interaction of the flow with topography, density-compensated submesoscale filaments (also called thermohaline filaments) induce the subduction of surface properties down into the water column (Fig. 3).
b. Feature definition: Spiciness
Our targeted “features” are thermohaline structures, subducting from below the mixed layer into the deep ocean (Fig. 3). Thermohaline filaments are a signature of strong horizontal and vertical mixing in the ocean (Klein et al. 1998). They may result from double-diffusion instabilities, as occurs in the Sargasso Sea (Ruddick and Gargett 2003), but they can also be the result of an efficient mechanism driven by mesoscale eddies to stir large-scale gradients into strong small-scale thermohaline fronts (Ferrari and Polzin 2005). Subduction regions are characterized by large kinetic energy and strain fields (Figs. 2 and 3).
It has been proposed that thermohaline filaments result from horizontal stirring of mesoscale temperature and salinity gradients, coupled with vertical advection by ageostrophic circulation (Smith and Ferrari 2009). This mechanism shows evidence of a direct cascade of mesoscale thermohaline variance to small (submeso) scales (Klein et al. 1998; Smith and Ferrari 2009).
Because these intrusions of near-surface waters are density compensated, they are found along constant density surfaces and therefore they become evident in spiciness π (Figs. 2g–i, 3c,f). The concept of a quantity related to isopycnal thermohaline variations and least correlated with the density field was first introduced by Stommel (1962), mathematically defined by Veronis (1972), and named “spice” (for warm and salty) by Munk (1981). Spiciness is useful in characterizing water masses and along-isopycnal mixing, and since the late 1980s it has been widely used to study the California Current System (Flament 2002). In this work, we use the lateral gradient of spiciness (/x; Fig. 3f) to detect features of interest.
Underwater gliders are suitable for the study of thermohaline filaments. Rudnick and Cole (2011) showed that as a result of the relatively slow nature of the glider, isobaric properties obtained from glider data are valid at wavelengths larger than 30 km, while isopycnal properties (like spiciness) are valid at the spatial resolution of the glider. This is because internal waves, which have an impact on isobaric surfaces, are intrinsically filtered out on isopycnal surfaces.
c. Underwater asset
Our ocean observing platform (or asset) is a Seaglider (Ogive profile; Kongsberg Underwater Technology, Inc.), an underwater autonomous buoyancy-driven vehicle capable of profiling to a maximum depth of 1000 m in a sawtooth pattern (V shape) (Eriksen et al. 2001). For the Monterey Bay experiment, the Seaglider carried a Sea-Bird SBE3 temperature sensor and SBE4 conductivity sensor (known collectively as the CT sail), a pressure sensor, and an Aanderaa 4330F oxygen optode. Following calibration, potential temperature, salinity, and oxygen concentrations are accurate to 0.01°C, 0.01, and 2 µmol kg−1, respectively. Sensor precision is 0.001°C and 0.0003 S m−1 for temperature and conductivity, respectively, combining to a salinity precision of approximately 0.001−1. Sampling occurred approximately every 5 s (0.5-m vertical resolution at typical vertical speeds of 0.1 m s). Compass calibration was performed before deployment.
The Seaglider data are processed using the University of East Anglia (United Kingdom) Seaglider Toolbox by B. Y. Queste (http://www.byqueste.com/toolbox.html), which includes tuning the hydrodynamic flight model following Frajka-Williams et al. (2011) and thermal lag corrections following the methods of Garau et al. (2011). Dive-average currents are calculated from the difference between the glider’s flight path found from GPS positions at the beginning and end of each dive, and the glider’s flight path as calculated from the Seaglider hydrodynamic model.
The methodology consists of a forecast numerical model that is used to predict the location of maximum lateral gradients of spiciness and an autonomous planner that commands the underwater asset. Shore-based feature detection software analyzes the model output and updates navigation commands on the autonomous planner to direct the underwater asset to the selected location. For this experiment adaptive sampling is considered once a day. The different components of the method are described next.
a. The Monterey Bay Coastal Ocean modeling system
The Monterey Bay Coastal Ocean modeling system is based on a nested data-assimilating version of the Regional Ocean Modeling System (ROMS) (Haidvogel et al. 2000; Shchepetkin and McWilliams 2005). The model configuration consists of an outermost domain covering the entire California coastal ocean from Ensenada, Mexico, to north of Crescent City, California, and extending approximately 1000 km offshore at a horizontal resolution of 3.3 km; an intermediate level domain with a resolution of 1.1 km covering the coast from Point Reyes to Point Conception to about 250 km offshore; and an innermost domain that encompasses the greater Monterey Bay region to about 75 km offshore with a horizontal resolution of approximately 300 m (Fig. 2). In the vertical there are 40 sigma levels used in each domain. The atmospheric forcing is derived from hourly output from operational forecasts performed with the NCEP 5-km North American model (NAM). Tidal forcing obtained from the TPXO.6 global barotropic tidal model (Egbert and Erofeeva 2002) is added through lateral boundary conditions. An essential component of this system is the new two-step multiscale three-dimensional variational data assimilation (MS-3DVAR) scheme used to generate the three-dimensional ocean state estimates, a generalization of the 3DVAR methodology of Li et al. (2008a,b) described in detail in Li et al. (2013, 2015a,b).
The system produces a single 2-day forecast of ocean conditions each day using a nowcast (or analysis of the ocean state) at 0300 UTC as initial condition. While the model also produces an additional nowcast every 6 h, the autonomous sampling described in this article uses only the initial nowcast and forecast. The nowcast-forecast system is run daily in near–real time. This means we aim to deliver the nowcasts and daily forecast within nine hours after the start time for which it is valid (i.e., the 0300 UTC nowcast and forecast are to be delivered by 1200 UTC). The system incorporates all available real-time streams of data, gathered in situ or remotely sensed (Table 1), including the California Institute of Technology (Caltech)’s Seaglider and the Naval Postgraduate School (NPS)’s Spray glider data gathered as part of the KISS field experiment. The temperature, salinity, and depth-averaged velocities obtained from these two gliders were generally available during every ROMS daily window. More details on the modeling system, including the data assimilation methodology and a validation of the operational results obtained with the outermost (3.3 km) nest, can be found in Chao et al. (2018).
b. Feature detection
To detect features of interest, we leverage the concepts from computer vision and graph theory research communities to perform background model subtraction and automatically detect irregularly shaped objects (Stauffer and Grimson 1999; Li et al. 2002). The idea can be summarized into a three-stage process: (i) generate a background model of a given field; (ii) separate the information into two types, background and foreground features; and (iii) derive a score based on detected features.
The algorithm starts by building a background model of a given field P (in our case, P = /x), based on a probability distribution function of P and a critical threshold (Pc), that will be used to separate background “noise” from foreground features. For this contribution the algorithm uses a fixed value of Pc, above which all absolute values of P > Pc will be considered potential foreground features; the rest will be considered as background (Fig. 4a). Then the algorithm separates P into one of two states: 0 for background regions and 1 for foreground features. This results in a binary mask (or background model). After the background is established, the algorithm begins to scan the entire field using a traversal method for detecting connected components (Felzenszwalb and Huttenlocher 2004). Two points, A and B, are defined as connected if there exists a path of points such that , and ; and are neighbors. In other words, a connected component is defined as a region of points that is continuous and has no breaks in space. The algorithm identifies each connected component as a distinct object. This is of particular use to track the evolution of an object (variations of area, value, and location in time of a particular feature). The method also defines a minimum of k-connected points, denoting that at least k points must be connected, to be considered an object (feature).
c. Path planner
The Seaglider is remotely controlled by the path planner, which produces navigation directives for the glider. The execution of these directives results in the asset following a template path (transect). A template path is a series of waypoints in latitude–longitude–time space, where time is the expected time of arrival at the waypoint. For this deployment, straight transects were used as the template path.
The path planner determines the glider target waypoint given to the glider control software to follow a template path (transect). This is not a simple task, given that the glider has a low velocity (approximately 0.25 m s−1), while ocean currents can be easily higher (0.50 m s−1) and variable in space and time. The path planner operates by considering a range of control directions near the desired end location (Branch et al. 2016). Glider motion is simulated along a fixed heading, to a specified depth, using a fixed glide slope and speed that were defined by the operational setting and empirically validated from prior deployment data. At fixed time intervals, the interpolated ROMS velocity at the latitude, longitude, depth, and time of the glider is added to the glider’s own velocity. The path planner assesses the predicted end position for each of the range of control actions (e.g., glider target waypoints) and selects the glider target waypoint that results in the simulation-predicted glider position that is closest to the desired glider location.
We assume the modeled glider motion has a fixed control velocity, which, combined with the glide slope, predicts a vertical and horizontal velocity. These “control” vertical and horizontal velocities are then added to the predicted horizontal current velocities from ROMS to derive the predicted position. As ROMS is a grid-based model, we linearly interpolate in x, y, z, and t to estimate the currents that the glider will encounter.
The abovementioned motion model is of limited complexity—a more accurate model would account for the change in buoyancy at the top and bottom of the dive to infer the vertical force and therefore the vertical motion and lateral motion (Eriksen et al. 2001). In our several deployments of predictive control, anecdotally, model inaccuracy is a much greater source of error than the motion model inaccuracy, as measured by predicted versus actual depth-averaged current. However, a more rigorous evaluation of these factors and incorporation of a higher-fidelity model are both excellent areas for future work.
d. Decision-making: Feature detection planner
The goal of the feature detection (FD) planner is to optimize the glider’s sampling strategy. This is accomplished by 1) using feature detection techniques (section 3b) to detect regions of high-density gradients and 2) modifying the glider’s path using the path planner (section 3c) to direct the glider toward these regions.
In this experiment we aim to maximize the number of features encountered by the underwater asset (Seaglider) along a given template transect (red line in Figs. 5b,d). The template transect is 32 km long and, with an average glider velocity of 0.9 km h−1, it takes an average of 35.5 h to complete. Using a maximum diving depth of 600 m, this leads to ~17 dives per template transect. We give the FD planner a relatively restricted job: each morning, based on the 24-h model forecast, it should choose the best path for the glider, along the template transect. In other words, the FD planner will decide when, in the next 24 h, the glider should turn around.
The first step of the FD method is the extraction and scoring of features from the ROMS output. In the present work, the FD planner analyzes the model output using two-dimensional fields. The FD planner considers all possible turnaround points for the day’s glider projected path. Then, for each scenario, it extracts the ROMS forecast at the glider’s expected position and time. Each possible plan (whether the glider should turn around immediately, after one dive, two dives, etc.) is extracted from the ROMS forecast, by assuming that the glider follows the transect template perfectly and that there are no surface currents. Decision points are limited, as the glider can turn around only on the surface. We do not let the glider turn within two dives of the end of the transect.
Each possible scenario (transect) is scored based on the threshold method defined in section 3b. The FD code detects all features that meet the following criteria: (i) the feature has a (threshold) spiciness gradient larger than 7 10−6 kg m−4 and (ii) the feature is located in the vertical between 100 and 600 m, which is the maximum diving depth of the glider in this experiment. The total score of the transect is the sum of the gradient of spiciness of all features that meet these two conditions. For convenience, the total score is divided by 10. Note that the scoring could also be applied to three-dimensional fields.
Finally, the transect that scored the best is selected for execution by the glider and the commands are sent to the planner. On the next surfacing, the new plan is uploaded to the glider. As an example, Fig. 5 shows two (out of nine) possible scenarios for 25 October. On that day the feature detection planner decided that the best-case scenario (the one with the highest scores) corresponded to the glider turning after two dives (Figs. 5a,b).
4. Operational experience with ROMS
The feature detection planner ran once a day from 22 to 30 October 2016. A total of 120 V-shaped dives were collected, sampling from the surface to 600 m. The daily timeline was as follows. At 0830 local time (LT) the model output was uploaded to the server. At 0930 LT the feature detection code ran and sent appropriate commands to the planner events file.
Out of a total of nine days of the experiment, the decision-making software based on ROMS output decided to turn the glider around on eight separate occasions (“best turn” in Table 2). Many of the decisions were correctly executed, but not all. To evaluate whether we took the right directive (where to turn around), we compare the scores obtained from the executed turn to those from the best turn for each day of feature detection. We find that the executed scores were about ~90% of the best scores. The estimated efficiency loss was due to different contingencies that prevented the planner from executing the selected turn. Such contingencies included communication problems with the Seaglider and with the model output, and scoring issues that were solved during the course of the experiment.
To evaluate the benefits of feature detection we would, ideally, compare the transects sampled by a glider in feature detection mode versus a glider that would sample in a “traditional” mode (e.g., moving continuously back and forth along a transect). However, since we had only one glider available at the time of our experiment, we follow an alternative method. Using the model output, for each day we compare the score from the transect that scored best with the score from the transect without a turn (e.g., Fig. 5). The no-turn transect is equivalent to sampling without feature detection. The difference between these two scores is the “gain,” which ranges from 39% to 59% with two exceptions. On 25 and 30 October, the gain was 165% and 118%, respectively. Overall, feature detection improves the scoring by 49.5% (Table 2; Fig. 6).
5. Comparison of glider and model results
a. Feature detection phase
Next we compare the ROMS output to glider data obtained during the feature detection experiment. To construct each section, we extract the ROMS output at the exact glider location (Fig. 7a) and time. Glider data are optimally interpolated (OI) onto a regular grid before calculating the lateral gradients of spiciness. Each glider transect is divided into 50 grid points, which results in a horizontal grid spacing of about 300 m (i.e., similar to the ROMS output). The OI scheme uses three main parameters to smooth the data: a given number of grid points in the horizontal direction (d, light gray lines), a given distance in the vertical direction (p, m; dark gray lines), and the relative error (s, for which 0 < s < 1; blue). The parameters chosen for this study (bold black line) are , , and , which represent horizontal length scales of 1.5 km and vertical length scales of 75 m (given a vertical grid spacing of 5 m). Sensitivity tests (for p = 10, 30; d = 3, 7, 10; and s = 0.05, 0.2) show that the scoring is least sensitive to vertical smoothing and most sensitive to horizontal smoothing (Fig. 8).
Comparisons between observed and modeled fields reveal mean modeled-observation misfits of 0.024°C and 0.003 salinity units. Positive misfits indicate the model is overall warmer and saltier than in situ observations. Standard deviations of the misfits are 1.787°C in temperature and 0.165 salinity units. These represent root-mean-square errors of the predicted-minus-observed temperatures and salinities, normalized by the root-mean-square of the observations and expressed as a percentage, of less than 25% in temperature and 0.5% in salinity. The mean misfit in spiciness is −0.026 kg m−3 with a standard deviation of 0.152 kg m−3. For lateral gradients of spiciness, the mean misfit between model and observations is −2.5 10−6 kg m−4 with a standard deviation of −1.1 10−6 kg . Normalization of the standard deviation of the misfits by the standard deviation of the observations leads to ratios of 1.35 for temperature, 1.32 for salinity, 1.25 for spiciness, and 1.68 for the lateral gradient of spiciness. It is fair to admit that these are worst-case errors. Since there are often spatial and temporal shifts between a model and observations derived from the turbulent nature of the ocean, single-point model evaluations can lead to an unrealistically poor assessment of the numerical model.
The probability distribution function of the lateral gradient of spiciness shows robust agreement between ROMS output and glider data, indicating that the model successfully reproduces the variability observed by the glider (Figs. 4c,d). A more detailed description of the system as well as a comprehensive validation of the 1-km model output can be found in Chao et al. (2018).
At the beginning of the feature detection phase, a strong anticyclonic eddy developed, and the feature detection planner commanded the glider to stay at the offshore side of the transect, where the gradients were stronger (Fig. 9c). Model output and glider data show subducting thermohaline structures from the near surface to 400-m depth characterized by lateral gradients of spiciness on the order of 10−5 kg m−4 (Figs. 9a,b). The vertical (~50–200 m thick) and lateral (~5–20 km) scales of these structures are similar in both datasets. The glider data show larger variability in the vertical, whereas the intrusions appear more coherent in the vertical in the ROMS output. Modification of the OI parameters can produce a more coherent vertical structure in the glider data; however, this smooths over real vertical finescale structure in the ocean, which is not resolved in the ROMS model. Nevertheless, this difference does not qualitatively impact the ability of our algorithm to identify frontal regions from the ROMS output.
Toward the end of the experiment, the largest gradients were found on the onshore side of the transect (Fig. 10c). When the glider moves toward these shallower regions, thermohaline structures are found close to the seafloor, extending up to ~200 m above the bottom (Figs. 10a,b). However, on one occasion, the glider data showed a subsurface filament sinking to near the bottom (not shown) that the model did not accurately reproduce. At present, we attribute this episode to possible limitations of the model to reproduce finescale bottom boundary layer (BBL) processes. Accurate reproduction of BBL processes is challenging for model forecasts because they require fine vertical resolutions and a reliable bathymetry. Our ROMS configuration shows large differences (up to 200 m) between the bathymetry measured by the glider altimeter and the ROMS bathymetry, in particular in regions of steep bathymetry gradients, such as near the shelf break. Such discrepancies are due to the necessary bathymetric smoothing required by ROMS to run without using an impractically small time step (Slørdal 1997), but they will likely have a strong impact on the model uncertainty to forecast the exact location of topographically induced submesoscale features.
b. Comparison with traditional sampling phase
In an attempt to further address the benefits of our method using in situ data (in section 4 we used numerical output only), next we present an analysis of the same transect used for feature tracking, only this time sampled in “traditional mode.” We refer to traditional mode sampling when the glider is repeating a transect based on target waypoints, and not based on feature detection methods. The following dataset was obtained in September 2016.
Five repeats were performed along the template transect from 12 to 26 September 2016 (Fig. 11a), with an average of 72 h to complete each transect. To compare this to scores obtained during the feature detection phase, each transect was divided into two subtransects of ~18 km each before scoring (Table 4).
During this traditional sampling phase, the maximum gradients of spiciness were located near the shelf break (Fig. 12). The scores were consistently larger close to the shelf, both in ROMS and glider data (Fig. 11b). During traditional sampling the glider was not able to choose where to sample; it simply kept repeating the template transect. This situation contrasts with the feature detection phase, during which a strong anticyclonic eddy developed, and the feature detection planner commanded the glider to stay on the offshore side of the transect, where the gradients were stronger. The main conclusion from these analyses is that our method works at selecting areas with a larger spiciness gradient.
While model output and glider data show different absolute scorings on each transect, they have comparable scoring variations along each transect. This behavior is observed during both (feature detection and traditional sampling) phases (Figs. 7b, 11b). The key point here is that both the glider and the model find enhanced features in the same regions, rather than their absolute score. It is not surprising that the ROMS and glider scores are not exactly the same, because there are processes happening in the ocean that are not resolved by numerical models. Even subsampling in situ data at the same resolution of the model will not get exactly the same score values. Our goal is to identify regions with a high probability of submesoscale activity and to direct our underwater assets toward it in an autonomous fashion. Therefore, for our method to work, we require that the spatial distribution of the scores show similar variation.
We would like to note that the number of turns and dives per day should not be taken literally to evaluate the benefits of our method. We test the concept in Monterey Bay for practical reasons, and with limited operating time. The main novelty of our method is the use of the numerical output to forecast the glider path (not only to guide it). By analyzing all the possible virtual glider paths for a given day, and selecting the optimal choice, our method takes advantage of the model forecast to detect which region of the ocean is more relevant to sample (according to a set of criteria). This approach could aid in exploratory field programs, in deep ocean surveys when surfacing happens less frequently, and in remote areas that are not well known. The objective of the Monterey Bay experiment is to show that a numerical model that assimilates glider data and data from other observing platforms can be used in an autonomous fashion to optimize and control an underwater vehicle sampling strategy.
6. Previous surveys in Monterey Bay
A number of other approaches have been used in the Monterey Bay region to identify or to track frontal structures. In particular, long time series composed of cross-slope transects have been collected as part of the California Cooperative Oceanic Fisheries Investigations (CalCOFI) since the early 1950s to study the variability of the CCS (McClatchie 2014). CalCOFI now employs underwater gliders to increase the spatial and temporal resolution over ship-based surveys (McClatchie 2014; Rudnick 2016). This dataset has been instrumental in providing a climatology of the CCS that includes a mean field, an annual cycle, and the anomaly from the annual cycle for each hydrographic variable (Rudnick et al. 2017), as well as examining interannual variations of the circulation, for instance, in response to El Niño and La Niña events (Todd et al. 2011a; Jacox et al. 2016; Rudnick et al. 2017) and the 2014–15 warming anomaly (Zaba and Rudnick 2016).
Over shorter temporal and spatial scales, these long-term glider datasets have also provided statistical descriptions of mesoscale and submesoscale fronts, filaments, and eddies of the CCS. Key results include a strong seasonality in near-surface fronts and an interesting vertical structure in lateral mixing as evidenced by the maxima and minima in spice variance along density surfaces (Davis et al. 2008; Todd et al. 2011b, 2012; Powell and Ohman 2015). However, because of the focus on fixed transects in CalCOFI, these surveys differ from our adaptive approach, which has the potential to provide additional information on the evolution of active submesoscale regions at subdaily time scales.
We postulate that, in addition to the CalCOFI glider sections, an additional fleet of gliders operated using the autonomous decision-making method presented here could be used to predict the evolution of the eddies and fronts, and to navigate the gliders to obtain finer scale sampling of these features. With the glider capacity to optimize its navigation to rapidly reach the desired waypoint while maintaining an optimal distance between assets (Davis et al. 2009), underwater glider networks appear to be one of the best approaches to achieving the subsurface spatial resolution necessary for ocean research (Davis et al. 2002; Rudnick et al. 2004; Rudnick 2016). Now that the future of ocean surveying seems to lead toward a combined cluster of platforms (Rudnick 2016), an autonomously driven fleet of underwater gliders would contribute with finescale sampling in those regions of the ocean where submesoscale variability concentrates.
We also note that our feature detection method differs from other feature-tracking techniques, often used to track fronts (Cruz and Matos 2010, 2014; Zhang et al. 2010, 2012, 2016; Sun et al. 2016) or oil spills (Zhang et al. 2011). These feature-tracking efforts focus on following a feature (in most cases a surface front) in near–real time. Here, we do not aim to track the exact location of a front. Instead, we seek to predict regions where submesoscale activity is likely to be active. To this end, we use a numerical model, in forecast mode, to evaluate a number of future glider sampling scenarios, to score them based on a given set of parameters, and to make navigation decisions based on these scores. Future surveys that may implement many tens of gliders at one time would be more feasible with both increased piloting autonomy and guidance from regional data-assimilating models.
7. Summary and remarks
This work presents a new method for autonomously sampling submesoscale ocean structures using feature detection techniques. The method relies on a real-time data-assimilative ROMS forecast, a path planner driving the underwater vehicle, and off-board autonomous decision-making. Critically, this approach moves toward a fully autonomous framework in which all parts of the navigation are conducted without human intervention (Thompson et al. 2017).
We evaluate the feasibility of the method to track submesoscale thermohaline (density compensated) features using a Seaglider in Monterey Bay. The model successfully reproduces similar features to those observed in nature (vertical and lateral scales, number of features per transect). The method successfully directs the glider to regions of high-density gradients. Based on model evaluations, the method allows for a sampling gain of 50% of the magnitude of isopycnal tracer variability. Differences between model output and the glider data are attributed to the interpolation of the glider data and model uncertainties in the exact reproduction (in both location and time) of submesoscale features.
From a scientific point of view, feature detection has the benefits of being able to follow the desired ocean structure and to study its evolution in time and space. Our two-dimensional approach (template transect) is not optimal, since there may be missed opportunities to find more features if the glider goes off the designated transect. Our method could be expanded into a three-dimensional analysis of the model output and a multivehicle approach. Modulation of the maximum depths of the dives would increase the spatial resolution of sampling so that the features can be better sampled with higher confidence. Finally, we need to carefully consider that scoring-based decisions will likely put emphasis on repetitive sampling (back and forth) of a local feature. This will work for studies involving the rapid evolution of bottom-intensified flows, local studies on mixed layer depth variability, and the evolution of surface fronts (Zhang et al. 2016). However, when the scientific objective includes a larger source of spatial variability (e.g., a comprehensive study on the evolution of processes along the edges of a large mesoscale eddy), the experimental design requires a multivehicle approach and the inclusion of multiple weighting scores to isolate the object of study from coexisting processes with similar signatures.
From an autonomy point of view, our experiment can be viewed as a closed-loop feature prediction (via the ocean model) and adaptive control to follow the best path to observe ocean features predicted by the ocean model. We acknowledge that the span of autonomy is modest—specifically that the only variable under control of the autonomy software is the number of dives executed before the glider turns around, which is a discrete set of choice points, with only two alternatives at each choice, with the glider position considered as a one-dimensional position. This simplified version of a more general problem enables a more focused study of the search space and evaluation of the scores in simulation. However, we acknowledge that this experiment represents only a simple first step toward applying the same autonomy approach to less constrained choice problems, for example, enabling the vehicle to prosecute features in the full lateral (x, y) domain, enabling variability in the glide slope of the glider, varying the dive maximum depth, and many others.
Looking to the future, the underwater vehicle should be directed based on a combination of recently acquired vehicle data and ROMS predicted data. While the ROMS model can predict and take into account many inputs, prediction is challenging and working from actual data can avoid issues with model inaccuracies. In addition, we plan to formulate a coordinated response based on the assimilation of data acquired from multiple vehicles to further improve sampling. A future goal of the program is to perform a coordinated multivehicle search to track the evolution of oceanographic features.
We thank Hongchun (Carrie) Zhang for her help with the ROMS modeling; Giuliana Viglione, Zack Erickson, and Xiaozhou Ruan for their help deploying, trimming, and recovering the Caltech glider; and Patrice Klein for his useful comments. Color maps used in this contribution are from Thyng et al. (2016). This work is funded by the Keck Institute for Space Studies (generously supported by the W. M. Keck Foundation) through the project “Science-Driven Autonomous and Heterogeneous Robotic Networks: A Vision for Future Ocean Observations” (http://www.kiss.caltech.edu/new_website/techdev/seafloor/seafloor.html). This work was, in part, performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The authors are grateful to three anonymous reviewers, whose comments helped improve the manuscript.