Abstract

Ocean observing systems (OOSs) constituted by moorings and gliders are becoming relevant in oceanographic and climate studies. In these observing networks, the temporal variability is captured by mooring observations, while the spatial variability is obtained from gliders sampling in the surrounding area. The advent of this observing capability brings the need to find optimal procedures to sample a given ocean region with a glider in the presence of a neighboring mooring, in order to maximize the information content of the data collected by this observing network. Different criteria (e.g., A, G, or E optimality) commonly used in the geosciences to obtain an optimum design lead to different sampling strategies. The question of which criterion performs better for optimal design in the marine environment remains open. This work investigates optimal procedures to sample a given ocean region with a glider in the presence of a mooring. Specifically, observing systems simulation experiments (OSSEs) are carried out in the Ligurian Sea (western Mediterranean) in August 2010 to study the different sampling strategies. Three criteria, which respectively aim at minimizing the trace (A optimal), maximum diagonal value (G optimal), and maximum eigenvalue of the error covariance matrix (E optimal), are considered. The resulting temperature field estimations are evaluated against a control field at 50-, 100-, and 150-m depth. The results indicate that the most appropriate strategy for environmental characterization using gliders employs the A optimal criterion, minimizing the mean uncertainty over the study area.

1. Introduction

The marine environment is an extremely complex dynamical system. Its turbulent nature implies strong interactions between a wide range of spatiotemporal scales (Kuragano and Kamachi 2000). Additionally, marine chemical and biological processes are strongly affected by the physical variability (Robinson et al. 1999; Gargett and Marra 2002; Griffiths et al. 2005). The high spatiotemporal physical variability and the relevance of the interactions between the physical, chemical, and biological fields make the study of the ocean difficult: first, it requires simultaneous measurements of the physical, chemical, and biological parameters; and second, these measurements must be done with an adequate spatiotemporal resolution.

Ocean observations are mainly based on sensors and observing platforms (Dickey 2003). The spatiotemporal resolution of ocean observations relates to the characteristics of the observing platform. An observing platform must be capable of providing sustained synoptic observations, that is, to map ocean structures at adequate spatial resolution faster than significant changes occur. This requirement is not realized by many traditional platforms of oceanographic sampling because of their physical, economic, and/or operational limits. For this reason, ocean observations are transforming from platform-based capabilities to networks of sensor nodes (Curtin et al. 1993; Curtin and Bellingham 2009). It is envisioned that future ocean observations will be done by heterogeneous ocean observing networks involving static nodes (moored profilers and bottom-mounted systems), nodes with uncontrolled motion (drifter buoys and profiling floats), and nodes with controlled motion [ships, autonomous underwater vehicles (AUVs), gliders, and autonomous surface vehicles (ASVs)]. This in situ observing system will then be complemented by remote sensing platforms/sensors (including acoustic, aerial, or space based). Present ocean observatories, such as the U.S. Integrated Ocean Observing System (IOOS) and the Australian Integrated Marine Observing System (IMOS), among others, are heading toward this observational concept.

The advent of ocean observing systems involving heterogeneous observing platforms creates new scientific and technological demands (Malone and Cole 2000), such as the optimal allocation of observational resources to maximize the information content of the oceanographic data collected by the network. Exploiting the observing capabilities of the different nodes must be found by designing optimum sampling strategies to allow for an accurate representation of oceanographic processes. These sampling strategies could adapt to the evolution of the environment and consider possible limitations resulting from the motion of part of the platforms in the network.

In the past, designs of oceanographic sampling strategies generally relied on simple statistical estimates (Davies 1991), tests of the experimental design (Le Traon and Hernandez 1992) or, more often, subjective criteria. Bretherton et al. (1976) provided the first structured approach to optimize the sampling performance of a network of oceanographic sensors. Stimulated by the need to design a current meter array, they used objective mapping to determine the array configuration that produced a mapped field with minimum mean square expected error. However, the work lacked a defined methodology to search for an optimal solution in the configuration space. Instead, designs were selected based upon experience and intuition, and their performance was evaluated in terms of the corresponding objective functions. Barth and Wunsch (1990) and Barth (1992) incorporated nonlinear optimization methods (simulated annealing and genetic algorithms, respectively) to search for the best array configuration in different oceanographic sampling problems. Specifically, the former investigated the best static configuration of an array of acoustic tomographic sensors, while the latter considered the sampling design required for the time-dependent problem of determining the concentration of a tracer in a flow field. The problem of optimizing the sampling performance of a network of drifters was considered by Hernandez et al. (1995). Unlike other oceanographic platforms, the uncontrolled motion of the drifters does not allow a high-level planning of the sampling procedure. Instead, optimizing the cast strategy basically relied on choosing how many drifters to use and where to launch them. When utilizing genetic algorithms for optimization, the sampling strategy for 25 drifters yielded those deployment positions that provided homogeneous drifter coverage of a portion of the Azores region in space and time.

Optimizing cast strategies with gliders, the motion of which is controlled, has received more attention recently in the scientific and engineering communities. Gliders are underwater autonomous vehicles that are designed to sample vast areas of the interior ocean (Stommel 1989). Gliders make use of buoyancy changes, their hydrodynamic shape, and small fins to carry out zigzag motions between the surface and a predetermined depth with a net horizontal displacement. The nominal horizontal speed is approximately 2 km h−1. Coastal versions of gliders are limited to operate between 10- and 200-m depth. Long autonomy at sea is the main advantage of this platform.

The problem of optimizing sampling strategies to optimal paths for a glider sampling in a shallow and turbulent environment was considered by Galea (1999). Sampling with networks of gliders was studied in detail by Rudnick et al. (2004) and Leonard et al. (2007). Specifically, both studies explored the problem of sampling an oceanographic field with a network of gliders that is constrained to keep fleet formation, that is, maintaining the same distance between them. Most recent studies have been focused on determining optimal sampling strategies to minimize uncertainties in ocean model predictions. In these cases, sampling designs are guided by cost functions based on predicted errors when observations from the network are assimilated into the numerical model (Robinson and Glenn 1999; Lermusiaux et al. 2004; Ramp et al. 2009; Smith et al. 2010; Heaney et al. 2007). Predicted current fields are of special interest to infer the feasibility of sampling designs with gliders. Because of the slow glider motion, adverse current fields could prevent them from reaching all the waypoints in the design. Thus, finding an optimal sampling also implies solving reachability issues.

While techniques to determine optimal sampling designs for a specific platform or a network of those platforms are relatively mature, the study of optimum sampling designs with networks of heterogeneous platforms that differ from one another on speed, motion controllability, sampling rate, etc., are scarce in the literature. Alvarez et al. (2007) considered the problem of optimally merging networks of profiling floats and gliders. Specifically, a genetic algorithm was employed to find optimal glider trajectories, together with an unevenly distributed network of floats, to get the best-sampled field. The cost function to minimize was defined in terms of the mean formal error obtained from an optimum interpolation scheme. More recently, Leonard et al. (2010) examined the performance of a network with two types of gliders (four deep Spray gliders, diving up to 1500 m, and six coastal Slocum gliders, diving up to 200 m) that had integrated coordinated control of the fleet with adaptive sampling.

Observing system simulation experiments (OSSEs) provide a methodology to assess and design optimum sampling strategies in multiplatform ocean observing systems (OOS). The approach relies on the outputs of an ocean model simulation of the area monitored by the OOS, which are treated as “truth.” Virtual observations from different ocean observing platforms in the OOS are then simulated from the model run and analyzed in the same manner as real data. The performance of a given sampling strategy can be assessed by a metric depending on the difference between the original model output and the analyzed field resulting from the virtual observations. OSSEs have been used in oceanography to find optimal locations for moorings in the tropical and meridional Atlantic (Hackert et al. 1998; Hirschi et al. 2003, respectively) and in the tropical Indian Ocean (Oke and Schiller 2007; Ballabrera-Poy et al. 2007). Another example of the use of OSSEs for the evaluation of a multiplatform OOS can be found in Mourre et al. (2006), considering sea level measurements from both satellites and tide gauges.

Monitoring an ocean region with moorings and gliders constitutes one of the simplest configurations of an OOS. Despite its simplicity, this network structure is becoming relevant in oceanographic and climate studies. In these cases, the temporal variability is captured by mooring observations while the spatial variability is obtained from the gliders sampling in the surrounding area. Since 2006, gliders routinely complement mooring observations of the Cape Verde Observatory, off northwest Africa (Wallace 2009). Gliders and moorings were also essential components in the assimilation studies conducted in the Autonomous Ocean Sampling Network (AOSN-II; Shulman et al. 2005) and in the framework of the observational program of the Australian IMOS (Jones et al. 2012).

Building upon these past efforts, this article investigates optimal procedures to sample a given ocean region with a glider in the presence of a mooring. Optimality is not only investigated in terms of finding the best sampling trajectory of the glider, but also in terms of different optimal criteria commonly employed in optimum design theory. The article is organized as follows: section 2 defines the problem considered in a mathematical context and introduces the optimal criteria most commonly used in experimental design theory. In section 3, the numerical formalism developed to find optimal sampling designs with a glider and a mooring is proposed. In section 4 these procedures are used to compute and compare sampling designs with different optimal criteria in a restricted marine area in the Ligurian Sea (western Mediterranean Sea) where a permanent Eulerian observatory exits. The evaluation of the performance of the different optimal designs is done through OSSEs involving a glider and an Eulerian observatory. Finally, results are discussed in section 5.

2. Mathematical formalism

A gliderlike vehicle and an Eulerian observatory are assumed to sample a given ocean region. Temperature and salinity are the main physical variables sampled by gliders. This measurement capability, together with their slow motion, makes gliders well suited to monitor mesoscale variability. Measurements are assumed to be synoptic and thus, no time dependence will be considered in the analysis. Procedures could be extended to deal with nonsynoptic measurements, but the computational complexity and information requirements about the physical processes in a region would be substantially increased without adding significant insight to the problem considered. Thus, the problem addressed here is to plan glider trajectories that, when combined with mooring observations, provide the best representation of the mesoscale field in a synoptic time scale.

Oceanographically, the adequacy of a glider sampling strategy can be measured in terms of the dynamical information that can be extracted from the measured field. This is obtained by estimating, from the data gathered at arbitrarily locations, the best values at grid points of a regular grid. A common procedure in oceanography to estimate the values of scalar fields at unobserved points is the Gauss–Markov smoothing [GMS; also known as optimal interpolation or objective analysis (Stein 1999)]. The sampled field is interpreted as a weakly stationary or second-order stationary process defined by known background field and covariance function , where the angled brackets denote the ensemble average). The approach relies on the locations where data are collected and the a priori knowledge or in situ estimation of the covariance of the sampled field to provide the best linear estimation of the average and variance of the field at given unsampled locations. Specifically, the value of the field at an unobserved location is linked to a set of measurements using a discrete linear regression model. Appropriate weights are obtained from minimizing the variance of the estimate with respect to the regression coefficients. The GMS analyzed field and covariance are obtained by

 
formula
 
formula

where and are state vectors of analyzed and background fields at gridpoint locations, respectively, is the vector of observations, is the observation matrix, that is, a matrix operator that interpolates from the regular grid to measurement locations, Σ is the background covariance matrix of the grid points, and is the gain matrix defined by

 
formula

The superscript T indicates the transpose operation, and is the observation error covariance matrix. The latter will be assumed diagonal. A variation of the optimal interpolation procedure, called ensemble optimal interpolation (EnOI), has been proposed by Evensen (2003). In this approach, numerical simulations of the oceanographic conditions in the region of interest are employed to determine the background statistical representation of the sampled field [ and in Eqs. (1) and (2)]. In the present work the ergodic hypothesis is assumed, that is, the average of simulated oceanographic processes over a relatively long period of time and the average over the statistical ensemble are considered equivalent. Therefore, and are estimated here based on a long time integration of the model (30 days, in this study).

Ideally, the best experimental design would consist of visiting each grid point of the defined regular grid during a synoptic time scale. This sampling strategy can be very demanding in time and/or resources, and thus it is not feasible in practice. Alternatively, glider trajectories can be planned in such a way that the estimations of the field at the grid points optimize a certain metric. Specifically, in Eq. (2) is a measure of the reliability of field estimations when observations are available at certain locations. Thus, observations could be designed to make as small as possible in a certain sense. A measure of the “smallness” of a matrix is usually described by a scalar function of the matrix. An optimal experimental design can then be generically defined by

 
formula

where is a scalar function and Ω is the set of admissible experiments or set of observations compatible with logistic constraints (number of platforms, mission time, etc). One of the frequently used criteria is called D optimal design (Silvey 1980), where the determinant of the error matrix is taken as the scalar function to minimize. Other composite measures are also possible. A design is A optimal if it minimizes the trace of with respect to the location of observations. Thus, A optimal designs are targeted to reduce the overall level of uncertainty over a region. An E optimal design minimizes the maximum eigenvalue of the resulting covariance matrix. Eigenvectors of correspond to coherent spatial patterns of variability, whose associated variance is provided by their corresponding eigenvalues (the so-called explained variance). Thus, E optimal designs minimize the variance of the worst estimated spatial pattern of variability. In other words, this design aims at avoiding the presence of a highly predominant error mode. Finally, a G optimal design minimizes the maximum diagonal value in . This design prefers locations which are most likely to contain extreme values. Both D and G optimality are equivalent for linear models (Kiefer and Wolfowitz 1960). This equivalence is convenient when the covariance matrix has a large condition number because the G optimality criterion allows a trivial numerical resolution. Different designs may lead to different optimal paths for the same problem depending on the particular spatiotemporal variability of the region under consideration. In general, A optimal designs result in more globally uniform sampling strategies, while G optimal designs lump observations in areas of high uncertainties. The generalities of E optimal designs are less intuitive. Comparative examples of different optimal designs can be found in Barnes (1989) and Muller (2007). Finally, it is important to note that optimal designs imply the minimization of a scalar function of the estimated variance with respect to sampling locations. This differs from the minimization performed in GMS, which is done with respect to regression coefficients. Thus, a best linear regression model is computed for each tentative sampling strategy resulting from the different optimum sampling designs.

3. Numerical formalism

Three main numerical procedures are required to minimize Eq. (4). The first one is the numerical discretization of the problem into forms that are solvable using linear algebra methods. To do so, a two-dimensional finite-element mesh regular grid is considered. Specifically, the physical domain is subdivided into a finite number of area elements (quadrilaterals) such that the intersection of two elements is an edge, a corner, or empty. A corner of an element cannot lie on the edge of an adjacent element. The corners of the elements are called the nodes. A continuous scalar function is then approximated in terms of the nodal values and certain interpolation functions within each element,

 
formula

where is the value of the function at node or grid point i, and represents the interpolation functions expressed in a local coordinate system {r, s} (Dhatt and Touzot 1984). Piecewise polynomials of low-order or spline functions are usually employed as the interpolation functions. The specific mathematical expressions of these functions are not replicated here becuase they can be found in textbooks about finite elements (e.g., Dhatt and Touzot 1984; Zienkiewicz and Taylor 1995). Notice that Eq. (5) defines a continuous piecewise linear approximation of the continuous function . An advantage of this encoding is that the computational demand depends on the total number of nodes considered in the mesh and not on the size of the dataset. Finite-element techniques are very convenient when observations are gathered by gliders or AUVs because they permit the processing of a huge amount of data with limited computational effort (Alvarez 2011). Discrete expressions of the vector of state variables, background field, and covariance are then defined by their values at the grid nodes, , , and , where N is the total number of grid points in the mesh. Similarly, the value of the background field at observation locations can be estimated from Eq. (5), resulting in an expression for matrix in Eq. (3).

The automatic generation of glider trajectories compatible with prescribed mission constraints is a second numerical procedure required to solve Eq. (4). In this study, mission constraints are defined by the area of operation, the number of gliders (Ng) and their initial locations , the glider nominal speed (V), the total mission time (TM), the navigation time between waypoints (TW), the surfacing time (TS), and the current field (VC) vertically averaged up to a given depth in the water column determined by the diving depth of the glider. A trajectory for the jth glider starting at is defined through a sequence of waypoints and is made by straight-line segments connecting any two adjacent waypoints . The jth glider must cover in TM and each segment in TW. Surfacing occurs at time intervals given by TS, resulting in a total of surfacing locations per segment. At each surfacing location the gliders reposition, estimate the local current field, and correct their heading. The overall sampling mission is defined by , the set of feasible waypoints for all gliders in the network.

For optimization purposes, a sampling mission is encoded by commanded headings at waypoints . The set of waypoints corresponding to commanded headings is obtained from the following procedure: tentative waypoints for each jth glider are initially defined by the sequential relationship with . Real glider trajectories may differ from those defined by tentative waypoints because of the effect of the current field . A realistic trajectory of each jth glider navigation between two consecutive waypoints is then computed by considering a point-like model of the platform subjected to a velocity field resulting from adding the nominal speed in the heading direction and the local current field . A heading correction is applied at every time interval determined by the surfacing time parameter TS. Specifically, at each surfacing the jth glider estimates the current field by comparing the simulated surfacing position with that expected from dead reckoning, where , where and are the location and velocity derived from previous surfacing point, respectively. The current field is then estimated by

 
formula

and it is employed by the glider to find the appropriate heading to reach the tentative waypoint

 
formula

The surfacing position may differ from the tentative waypoint after the navigation time TW. In this case, is updated by this surfacing location and the process is repeated for the next navigation segment . Thus, final waypoints defining the mission are determined from the surfacing locations at each time interval defined by the navigation time parameter. The procedure ensures that the final set of waypoints are timely achievable by the platforms. In the present study, glider observations are located every half kilometer along the simulated glider trajectory. This would correspond to the horizontal sampling resolution of a coastal Slocum glider.

Discretization and simulation of glider navigation allow the computation of Eq. (3) for an arbitrary mission encoded by a set of headings. Equation (4) must then be solved for a selected optimal criterion. This requires solving a nonlinear optimization problem, with the nonlinearity being introduced by the fact that the a posterior covariance is nonlinear with respect to the observation matrix , which depends on sampling locations, through Eqs. (2) and (3). Such problems are difficult to solve when the cost function is too complicated or not differentiable. A best design may be found by a trial-and-error method, that is, by generating numerous feasible designs and evaluating their performance in terms of the optimal criterion selected. However, this approach does not guarantee that the global minimum will be located. Soft-computing optimization methods are an alternative to the trial-and-error method. Soft computing embeds those techniques that exploit tolerance for precision and uncertainty to provide approximate but timely solutions to computationally hard problems. In essence, these approaches are built on deterministic and/or stochastic searching strategies in the solution space. In general, their computational complexity only increases linearly with the dimension of the solution space. Their drawback is that convergence to the optimal solution is not guaranteed in a finite number of iterations, so that one may end up with a suboptimal solution. Pattern search optimization (PSO; Hooke and Jeeves 1961) is a soft-computing approach that attempts to get closer to the optimal solution, searching in a mesh built around a seed point with a given pattern. The latter is defined as a set of vectors that determine the searching directions. If a grid point improves the value of the cost function at the current point, the new point is taken as a temporary solution and a new searching mesh is generated around it. Conversely, the searching mesh is stretched along pattern directions if no grid point improves the current value of the cost function. The operation is repeated until certain stopping criteria are met. PSO was found to provide better performance than genetic algorithms and simulated annealing when searching for optimum sampling designs of a glider fleet (Alvarez and Mourre 2011). Thus, this approach has been selected here to solve the nonlinear optimization problem. In this study, the searching pattern is the “positive basis 2N” [Coope and Price (2002); with N and being the number of variables in the cost function and the ith unit basis in , respectively) and the mesh contraction and expansion factors are fixed to 0.5 and 2, respectively.

4. Optimum mission planning for a glider–mooring network

a. Method

The mathematical and numerical procedures described above have been employed to compute optimum mission design for a glider operating in a marine area where a mooring is also present. Specifically, the objective is to find an optimum glider trajectory to get near-optimal temperature estimations in the first 200 m of the water column. This depth was determined on the basis of the maximum diving depth of a Slocum coastal glider. The region selected was a rectangle of approximately 60 km × 60 km in the Ligurian Sea centered on the Ocean Data Acquisition System (ODAS) Italia 1-W1M3A Eulerian observatory (Fig. 1). The observatory, located at 43.79°N, 9.16°E, is instrumented by a buoy station and a close mooring, including acoustic Doppler current profilers (ADCPs) and conductivity–temperature–depth (CTD) sensors. Hypothetical ODAS mooring observations of the local temperature field up to 200-m depth are considered in this study.

Fig. 1.

Geographical location of the ODAS Italia 1-WIM3A Eulerian observatory (black circle), area of glider operations (gray square), and circulation pattern in the Ligurian Sea (arrows).

Fig. 1.

Geographical location of the ODAS Italia 1-WIM3A Eulerian observatory (black circle), area of glider operations (gray square), and circulation pattern in the Ligurian Sea (arrows).

Oceanographically, the ODAS mooring is located in one of the most dynamically active and productive areas within the entire Mediterranean. A permanent cyclonic circulation with noticeable seasonal and interannual variability dominates the area (Millot 1987). The cyclonic circulation affects the hydrological structure of the basin, giving rise to three distinct zones: the coastal peripheral, the frontal, and the central zones. The frontal zone divides the warmer and lighter coastal waters from the colder and denser waters of the central zone (Sournia et al. 1990). Although this structure is almost permanent, it shows an important seasonal and interannual variability and intense mesoscale activity related to the strong frontal system (Astraldi and Gasparini 1992; Vignudelli et al. 2000). A main current system, the so-called Northern Current (NC), flows along the Ligurian and French coasts (Fig. 1). The synoptic time scale in the region estimated from the model run is approximately 4–5 days.

Optimum mission designs for a glider were computed for the described area for the period from 20 to 24 August 2010, using A, G, and E optimal criteria. A 30-day integration during August 2010 of the operational Navy Coastal Ocean Model (NCOM; Martin 2000) is considered in this study. NCOM was coupled to the Navy Coupled Ocean Data Assimilation (NCODA) system (Cummings 2005), and was run at the U.S. Naval Research Laboratory (NRL) Stennis Space Center (SSC). Satellite sea surface temperature and height, as well as in situ data from CTD stations and gliders, were assimilated in the model. The model configuration has a 1.8-km horizontal resolution and a vertical resolution varying from 2 m at the surface to 10 m at 200-m depth. The time-evolving thermal field resulting from the 30-day simulation is considered as the “truth” from which glider and mooring data are simulated, and against which field estimations are evaluated. This experimental setup allows for ideal conditions for several experiments in a perfect model environment, and thus for concentration on the capability of field estimation using different optimal criteria rather than on intrinsic model errors.

An average temperature background field and associated covariance matrix were generated at control depths of 50, 100, and 150 m from the time series of daily average NCOM analysis. In addition, a vertically integrated covariance matrix was also calculated. The single-depth fields are used in the evaluation of the reconstructed temperature estimate, while the integrated covariance is involved in the glider mission optimization. The glider mission indeed intends to optimize temperature estimations in the first 200 m of the water column and not for a specific depth. However, the final field estimations are expected to be more accurate if single-depth background fields and covariances are employed in Eqs. (1), (2), and (3) instead of vertically integrated values. The observation error in (3) was fixed to 0.1°C. This value was estimated from the model on the basis of the variability found in the thermal field at the synoptic time scale. Thus, the spatiotemporal variability of the thermal field at scales faster than the synoptic time scale is incorporated into the formalism as a representation error of the platforms measurements. This parameterization, which is based on the spatiotemporal variability described by the numerical model, is acceptable in the numerical OSSE framework. However, this error should be increased in real situations to account for the wider range of ocean variability projected onto real glider observations (Rudnick and Cole 2011). The model numerical grid determines the grid points where the temperature field is estimated from mooring and glider observations. The 4-day sampling period was split into two mission phases of 2 days each. Phases were sequentially linked, that is, the final waypoint of the first phase is used as starting point for the second phase, and the updated covariance matrix is passed to the next cycle. This splitting technique allows for incrementing the resolution of the mission without increasing the dimensionality of the solution space, thus facilitating the optimization process. Respectively, TW and TS were fixed to 12- and 3-h periods, with a glider nominal speed of 0.35 m s−1.

Generally, results will be reported for the control depth of 50 m. Figures 2a,b display the mean temperature field and standard deviation at 50-m depth obtained from the monthly model run. The region is characterized by a 2°C thermal gradient from coastal to open waters. The variability field, used in our approach to approximate the model uncertainty, exhibits a marked frontal structure, with larger values in the southern part of the domain.

Fig. 2.

(a) Mean temperature field and (b) associated standard deviation at 50-m depth computed from the operational NCOM model during August 2010.

Fig. 2.

(a) Mean temperature field and (b) associated standard deviation at 50-m depth computed from the operational NCOM model during August 2010.

b. Results

The ocean temperature at 50-m depth and the current field integrated in the first 200 m are respectively displayed in Figs. 3a,b for the 20–24 August period. The coastal to open sea thermal gradient observed in the mean monthly field is present, with enhanced spatial structure derived from frontal instabilities. Current field is characterized by two main northward streams: while the first one is associated with the continental slope, the second corresponds to the cyclonic circulation observed in the center of the Ligurian Sea. Maximum vertically integrated speeds are 0.15 m s−1.

Fig. 3.

(a) Temperature field at 50-m depth and (b) current field integrated in the first 200 m of the water column derived from the operational NCOM model for 20–24 Aug 2010.

Fig. 3.

(a) Temperature field at 50-m depth and (b) current field integrated in the first 200 m of the water column derived from the operational NCOM model for 20–24 Aug 2010.

Figures 4a,b display, together with the glider trajectory, the estimated temperature and absolute error fields at 50-m depth obtained from simulated glider and mooring observations for the A optimal design. The glider surrounds the mooring to the south, before heading northward along the western boundary of the domain. The estimated temperature field nicely reproduces the model temperature field for that period. Absolute errors, defined as the absolute differences between the estimated and true fields, are lower than 0.2°C.

Fig. 4.

(a) The A optimal temperature field at 50-m depth and (b) associated error, (c) G optimal temperature field at 50-m depth and (d) associated error, and (e) E optimal temperature field at 50-m depth and (f) associated error for 20–24 Aug 2010. The glider trajectory and mooring location are displayed by the black line and point, respectively.

Fig. 4.

(a) The A optimal temperature field at 50-m depth and (b) associated error, (c) G optimal temperature field at 50-m depth and (d) associated error, and (e) E optimal temperature field at 50-m depth and (f) associated error for 20–24 Aug 2010. The glider trajectory and mooring location are displayed by the black line and point, respectively.

Results from the G optimal design are displayed in Figs. 4c,d. Similar to the A optimal design, the glider heads southward to turn around the mooring, until it approaches the southern boundary. Notice that the largest initial uncertainty is found in the southwest corner of the region. Therefore, the G optimal design drives the glider to visit those points. The glider is then directed northward, and the structure of the absolute error field is similar to the one obtained from the A optimal design. However, the magnitude of the error is larger in the G design scenario (with absolute values exceeding 0.2°C).

Figures 4e,f show the results obtained for an E optimal design. The resulting trajectory is less intuitive than in the previous cases because the E optimality attempts to minimize the magnitude of the main pattern of error variability, which may change with the selected observations. Performance of this optimal criterion is poorer than the previous designs. The maximum absolute error amplitudes are around 0.5°C. Interestingly, all designs share a similar spatial structure of the absolute error field, characterized by larger values in the southern part of the domain. Further analysis has demonstrated that this spatial structure originates from the intrusion in the domain of cold water masses from the open sea between 22 and 24 August, the period during which the glider is sampling the western and northern parts of the domain. As a consequence, the temperature field is overestimated south of the ODAS mooring during the second part of the 4-day sampling period. In this case, the time scales of variability are shorter than the synoptic scale and responsible for the observed spatial error structure.

A similar performance analysis has been carried out for the remaining control depths. Table 1 evaluates the estimated fields in terms of the root-mean-square error (RMSE). RMSE is commonly used as a sample statistic to assess the match between an estimation and observations. Results indicate that A optimal design leads to a better performance than G and E optimal sampling strategies. Notice that the differences between the sampling scenarios are reduced at 100 and 150 m because of the weaker temperature variability at these depths.

Table 1.

RMSE of the temperature estimations at 50, 100, and 150 m obtained from different optimal designs.

RMSE of the temperature estimations at 50, 100, and 150 m obtained from different optimal designs.
RMSE of the temperature estimations at 50, 100, and 150 m obtained from different optimal designs.

The chi-square test has been applied to the distribution of the residuals over the whole area to provide insights into the consistency between the statistical assumptions of the method and the resulting field estimates. Under the present statistical assumptions, the sum of the square of the residuals normalized by observation variance should follow a chi-square statistical distribution (Aster et al. 2005). The p value of the chi-square test indicates how independent and normally distributed the normalized residuals are. The p values very close to either zero (e.g., 10−12) or one indicate the failure of one or more assumptions (Aster et al. 2005). Respectively, p values for the A, G, and E optimal designs at 50-m depth are 0.008, 3 × 10−14, and 0, respectively. These results indicate a certain consistency between estimation and statistical assumptions for the A optimal design. This consistency is lost for the sampling configuration described in the G and E optimal designs. Notice that all of the designs share the same statistical assumptions and estimation approaches, and they only differ on the location of the observations. This loss of consistency highlights the presence in the modeling area of outlier values of field–observation differences when G and E designs are considered. Simulations were also done considering glider deployments at corner and middle–edge points of the rectangular boundary. The A optimal designs provided the best field estimations in 62% of the cases. In the remaining cases, the A optimal designs were always second ranked with small differences in RMSE with respect to the first one. However, the case previously reported is the only one where statistical consistency was achieved by any of the optimal designs.

5. Discussion

Allocating and complementing observational resources to maximize the information content of the collected data has gained interest in different fields of geophysical sciences, including oceanography. The objective is to get the best field estimations on a regular grid from observations gathered at sparse locations. Different criteria have been proposed in the literature to measure the optimality of estimations. The A optimality defines the best estimation as the one that minimizes the associated average uncertainty; G optimality measures the goodness of the estimation by the magnitude of the maximum value of this uncertainty; finally, E optimality focuses on minimizing the effect of the main spatial pattern of variability of the error. The question of which of these optimal criteria is more suited to estimate oceanographic fields has been investigated in the present study in the context of an heterogeneous ocean observing network formed by a glider and a mooring in a restricted 60 km × 60 km marine area. Our results suggest that, in the area considered, A optimality is more suited than G and E optimal designs for oceanographic estimations with a glider and mooring. In particular, the G optimal design is more sensitive to pointlike anomalies in the uncertainty field. This may bias the sampling strategy toward remote locations with high uncertainty, thus degrading the overall estimation. The E optimal design is discarded on the basis of our results.

The methodology developed to conduct this research involved procedures to (i) get field estimations at certain points of interest from observations gathered at arbitrary locations, (ii) generate feasible glider trajectories compatible with operational constraints, and (iii) minimize a cost function defined by a selected optimal criterion. Various techniques could be employed for each procedure. Here, the GMS or optimal interpolation has been chosen as the procedure used for field estimation. The approach relies on a Gaussian statistical description of the environment where a covariance function encodes the spatiotemporal variability in the area. Estimations are then obtained from the statistical model assuming synoptic observations; that is, it is assumed that measurements are done before the environment changes significantly. The synoptic assumption is of common use in observational oceanography because it simplifies field estimation methodologies by limiting the required amount of information about the physical processes in the region. This is a desirable property for the optimization of sampling designs, where field estimations need to be evaluated for a large number of tentative sampling configurations. More complex approaches could imply a significant computer demand and make optimization unfeasible. In reality, the marine environment is permanently evolving and the synopticity assumption may induce errors in the field estimate. In this study, the spatiotemporal variability within the synoptic scale has been modeled as a representation error of the observations, as a compromise to reduce the impact of nonsynopticity without increasing the complexity of the statistical description. This description is supported by the fact that the resulting representation error is significantly smaller (by a factor 3) than the total uncertainty of the field represented in the background covariance matrix. Notice that this does not hold for situations where the sampling period is longer than the synoptic time scale because the representation error of a measurement would converge to the total variance of the field. Thus, observations would not significantly improve the “historical” information encoded in the background covariance matrix. More sophisticated and time-dependent descriptions of the environmental variability would then be required to find optimum sampling strategies of the proposed heterogeneous ocean observing network.

A code was developed to generate glider missions compatible with operational constraints. An important aspect concerns the oceanic current field. Gliders strongly interact with currents because of their slow speed. Waypoints solely determined on the basis of the glider nominal speed and navigation time could not be realistic because of the existence of adverse currents. To anticipate the impact of the currents, a mission planner should consider current field estimations. The accuracy of the current estimations considered by the mission planner may significantly impact the field performance of the different sampling strategies. The quantification of this impact needs to be investigated, and it still remains an open question.

A pattern search optimization (PSO) technique was selected to minimize the cost function. The selection was based on authors’ experience during field experiments. Specifically, PSO showed a better performance than genetic algorithms (GA) and simulated annealing (SA) for a relatively short computation time (600 s). This does not indicate a superiority of PSO versus GA and SA, but instead illustrated a faster convergence toward a near-optimum solution. Notice that the stochastic nature of GA and SA makes these algorithms more robust to local minima while it slows down their convergence. It is expected that GA and SA would get better performance than PSO for long computation periods. However, getting a near-optimal solution in a relatively short time period is preferable in operational applications.

The performance assessment of the sampling strategies resulting from the different optimality criteria raises a particular issue concerning the criterion used to measure this performance. The RMSE is a well-established criterion used to evaluate the realism of a field estimation. The fact that the A optimality seeks sampling designs that minimize the ensemble average (or the time average if ergodicity is assumed) of the RMSE raises the question of whether the choice of the RMSE criterion biased the field evaluations in this study by favoring the A optimal design. Concerning this question, let us first mention that the minimization of the ensemble average of a given variable does not imply the minimization of this variable for a single realization of the statistics, which has been considered here to evaluate the field estimate. Second, in addition to its numerical value, the sum of the squares of the residuals also provides useful statistical information about the quality of the estimates. The chi-square test, which measures the consistency of the statistical assumptions and mathematical models used in the estimation, has only revealed statistical consistency in the case of the A optimal design. Summarizing, even if our results are inevitably linked to the choice of the RMSE as the measure of the performance, and to the regional ocean dynamics, they indicate that the A optimal design provided both the most real (in an RMSE sense) and the most statistically consistent estimation of the thermal field in the study area.

Finally, the estimation error field was shown to be impacted by the violation of the synoptic assumption in this study, independent of the optimization criterion applied. This aspect is related to the proximity of a frontal structure, which generates a significant temporal variability at time scales shorter than the assumed synoptic scale of 4 days. As a consequence, the dimension of the area to be covered by a single glider during this time scale should be reduced. Alternatively, networking several glider platforms would constitute a suitable procedure to overcome this problem and get fully synoptic descriptions of spatially extended marine areas. Potential procedures to assess the minimum number of gliders required to synoptically monitor a given area are presently being investigated.

Acknowledgments

This work has been partially supported by the project EU-FP7-SST-234096-Argomarine and the NURC Program of Work. Authors want to thank Emanuel Coelho and Germana Peggion from NRL Stennis Space Center.

REFERENCES

REFERENCES
Alvarez
,
A.
,
2011
:
Volumetric reconstruction of oceanographic fields estimated from remote sensing and in situ observations from AUVs of opportunity
.
IEEE J. Oceanic Eng.
,
36
,
12
24
.
Alvarez
,
A.
, and
B.
Mourre
,
2011
:
Optimum and adaptive mission planning of gliders. Proc. Fifth European Glider Observatories Meeting, Canary Island, Spain, European/Everyone’s Gliding Observatories Network. [Available on http://ego2011.eu/doc/PDF-EGO/Friday%2018th/EGO%202011%20-%20Alberto%20Alvarez%201.pdf.]
Alvarez
,
A.
,
B.
Garau
, and
A.
Caiti
,
2007
:
Combining networks of drifting profiling floats and gliders for adaptive sampling of the ocean. Proc. 2007 IEEE Int. Conf. on Robotics and Automation, Rome, Italy, IEEE, 157–162
.
Aster
,
R. C.
,
B.
Borchers
, and
C. H.
Thurber
,
2005
:
Parameter Estimation and Inverse Problems. Elsevier Academic Press, 291 pp
.
Astraldi
,
M.
, and
G. P.
Gasparini
,
1992
:
The seasonal characteristics of the circulation in the north Mediterranean basin and their relationship with the atmospheric-climatic conditions
.
J. Geophys. Res.
,
97
(
C6
),
9531
9540
.
Ballabrera-Poy
,
J.
,
E.
Hackert
,
R.
Murtugudde
, and
A. J.
Busalacchi
,
2007
:
An observing system simulation experiment for an optimal moored instrument array in the tropical Indian Ocean
.
J. Climate
,
20
,
3284
3299
.
Barnes
,
R. J.
,
1989
:
Sampling design for geologic site characterization. Geostatistics: Proceedings of the Third International Geostatistics Congress, Vol. 2, M. Amstrong, Ed., Kluwer Academic Publishers, 809–822
.
Barth
,
N.
,
1992
:
Oceanographic experiment design II: Genetic algorithms
.
J. Atmos. Oceanic Technol.
,
9
,
434
443
.
Barth
,
N.
, and
C.
Wunsch
,
1990
:
Oceanographic experiment design by simulated annealing
.
J. Atmos. Oceanic Technol.
,
20
,
1249
1263
.
Bretherton
,
F. P.
,
R. B.
Davies
, and
C. E.
Fandry
,
1976
:
A technique for objective analysis and design of oceanographic experiments applied to MODE-73
.
Deep-Sea Res.
,
23
,
582
599
.
Coope
,
I. D.
, and
C. J.
Price
,
2002
:
Positive bases in numerical optimization
.
Comput. Optim. Appl.
,
21
,
169
175
.
Cummings
,
J.
,
2005
:
Operational multivariate ocean data assimilation
.
Quart. J. Roy. Meteor. Soc.
,
131
,
3583
3604
.
Curtin
,
T.
, and
J.-G.
Bellingham
,
2009
:
Progress towards autonomous ocean sampling networks
.
Deep-Sea Res. II
,
56
,
62
67
.
Curtin
,
T.
,
J.-G.
Bellingham
,
J.
Catipovic
, and
D.
Webb
,
1993
:
Autonomous oceanographic sampling networks
.
Oceanography
,
6
(
3
),
86
94
.
Davies
,
R. E.
,
1991
:
Observing the general circulation with floats
.
Deep-Sea Res.
,
38
,
531
571
.
Dhatt
,
G.
, and
G.
Touzot
,
1984
:
The Finite Element Displayed. John Wiley & Sons, 509 pp
.
Dickey
,
T.-D.
,
2003
:
Emerging ocean observations for interdisciplinary data assimilation systems
.
J. Mar. Syst.
,
40–41C
,
5
48
.
Evensen
,
G.
,
2003
:
The ensemble Kalman filter: Theoretical formulation and practical implementation
.
Ocean Dyn.
,
53
,
343
367
.
Galea
,
A. M.
,
1999
:
Various methods for obtaining the optimal path for a glider vehicle in shallow water and high currents. Proc. Int. Conf. Unmanned Untethered Submersible Technology, Durham, NH, Autonomous Undersea System Institute, 150–161
.
Gargett
,
A. E.
, and
J.
Marra
,
2002
:
Effects of upper ocean physical processes turbulence, advection, and air-sea interaction on oceanic primary production. The Sea, A. Robinson, B. Rothschild, and J. McCarthy, Eds., Biological-Physical Interactions in the Sea, Vol. 12, John Wiley & Sons, 19–49
.
Griffiths
,
G.
,
S.
Fielding
, and
H. S. J.
Roe
,
2005
:
Biological–physical–acoustical interactions, The Sea, A. Robinson, B. Rothschild, and J. McCarthy, Eds., Biological-Physical Interactions in the Sea, Vol. 12, John Wiley & Sons, 441–474
.
Hackert
,
E. C.
,
R. N.
Miller
, and
A. J.
Busalacchi
,
1998
:
An optimized design for a moored instrument array in the tropical Atlantic Ocean
.
J. Geophys. Res.
,
103
,
7491
7509
.
Heaney
,
K. D.
,
G.
Gawarkiewicz
,
T. F.
Duda
, and
P. F. J.
Lermusiaux
,
2007
:
Non-linear optimization of autonomous undersea vehicle sampling strategies for oceanographic data-assimilation
.
J. Field Robot.
,
24
(
Special Issue
),
437
448
.
Hernandez
,
F.
,
P. Y.
Le Traon
, and
N. H.
Barth
,
1995
:
Optimizing a drifter cast strategy with a genetic algorithm
.
J. Atmos. Oceanic Technol.
,
12
,
330
345
.
Hirschi
,
J.
,
J.
Baehr
,
J.
Marotzke
,
J.
Stark
,
S.
Cunningham
, and
J. O.
Beismann
,
2003
:
A monitoring design for the Atlantic meridional overturning circulation
.
Geophys. Res. Lett.
,
30
,
1413
,
doi:10.1029/2002GL016776
.
Hooke
,
R.
, and
T. A.
Jeeves
,
1961
:
Direct search solution of numerical and statistical problems
.
J. Assoc. Comp. Mach.
,
8
,
212
229
.
Jones
,
E. M.
,
P. R.
Oke
,
F.
Rizwi
, and
L.
Murray
,
2012
:
Assimilation of glider and mooring data into a coastal ocean model
.
Ocean Modell.
,
47
,
1
13
,
doi:10.1016/j.ocemod.2011.12.009
.
Kiefer
,
J.
, and
J.
Wofowitz
,
1960
:
The equivalence of two extremum problems
.
Can. J. Math.
,
12
,
363
366
.
Kuragano
,
T.
, and
M.
Kamachi
,
2000
:
Global statistical space-time scales of oceanic variability estimated from the TOPEX/POSEIDON altimeter data
.
J. Geophys. Res.
,
105
,
955
974
.
Leonard
,
N. E.
,
D.
Paley
,
F.
Lekien
,
R.
Sepulchre
,
D. M.
Frantantoni
, and
R.
Davis
,
2007
:
Collective motion, sensor networks and ocean sampling
.
Proc. IEEE
,
95
,
48
74
.
Leonard
,
N. E.
,
D.
Paley
,
R. E.
Davis
,
D. M.
Fratantoni
,
F.
Lekien
, and
F.
Zhang
,
2010
:
Coordinated control of an underwater glider fleet in an adaptive ocean sampling field experiment in Monterey Bay
.
J. Field Robot.
,
27
,
718
740
.
Lermusiaux
,
P. F. J.
,
C.
Evangelinos
,
R.
Tian
,
P. J.
Haley
,
J. J.
McCarthy
,
N. M.
Patrikalakis
,
A. R.
Robinson
, and
H.
Schmidt
,
2004
:
Adaptive coupled physical and biogeochemical ocean predictions: A conceptual basis. Computational Science— ICCS 2004: Proceedings of the Fourth International Conference, F. Darema, Ed., Lecture Notes in Computer Science, Vol. 3038, Springer, 685–692
.
Le Traon
,
P. Y.
, and
F.
Hernandez
,
1992
:
Mapping the oceanic mesoscale circulation: Validation of satellite altimetry using surface drifters
.
J. Atmos. Oceanic Technol.
,
9
,
687
698
.
Malone
,
T. C.
, and
M.
Cole
,
2000
:
Toward a global scale coastal ocean observing system
.
Oceanography
,
13
,
7
11
.
Martin
,
P. J.
,
2000
:
Description of the Navy Coastal Ocean Model Version 1.0. Naval Research Laboratory Rep. NRL/FR/7322-00-9962, 42 pp
.
Millot
,
C.
,
1987
:
Circulation in the Western Mediterranean Sea
.
Oceanol. Acta
,
10
,
143
149
.
Mourre
,
B.
,
P.
De Mey
,
Y.
Menard
,
F.
Lyard
, and
C.
Le Provost
,
2006
:
Relative performance of future altimeter systems and tide gauges in constraining a model of North Sea high-frequency barotropic dynamics
.
Ocean Dyn.
,
56
(
5–6
),
473
486
.
Muller
,
W. G.
,
2007
:
Collecting Spatial Data. Springer, 242 pp
.
Oke
,
P. R.
, and
A.
Schiller
,
2007
:
A model-based assessment and design of a tropical Indian Ocean mooring array
.
J. Climate
,
20
,
3269
3283
.
Ramp
,
S.
, and
Coauthors
,
2009
:
Preparing to predict: The Second Autonomous Ocean Sampling Network (AOSN-II) experiment in Monterey Bay
.
Deep-Sea Res. II
,
56
,
68
86
.
Robinson
,
A. R.
, and
S. M.
Glenn
,
1999
:
Adaptive sampling for ocean forecasting
.
Nav. Res. Rev.
,
51
,
28
38
.
Robinson
,
A. R.
,
J. J.
McCarthy
, and
B. J.
Rothschild
,
1999
:
Interdisciplinary ocean science is evolving and a systems approach is essential
.
J. Mar. Syst.
,
22
,
231
239
.
Rudnick
,
D. L.
, and
S. T.
Cole
,
2011
:
On sampling the ocean using underwater gliders
.
J. Geophys. Res.
,
116
,
C08010
,
doi:10.1029/2010JC006849
.
Rudnick
,
D. L.
,
R. E.
Davis
,
C. C.
Eriksen
,
D. M.
Frantantoni
, and
M. J.
Perry
,
2004
:
Underwater gliders for ocean research
.
Mar. Technol. Soc. J.
,
38
,
48
59
.
Shulman
,
I.
,
D. J.
McGillicuddy
,
M. A.
Moline
,
S. H. D.
Haddock
,
J. C.
Kindle
,
D.
Nechaev
, and
M. W.
Phelps
,
2005
:
Bioluminescence intensity modeling and sampling strategy optimization
.
J. Atmos. Oceanic Eng.
,
1267
1281
.
Silvey
,
S. D.
,
1980
:
Optimal Design—An Introduction to the Theory for Parameter Estimation. Chapman and Hall, 86 pp
.
Smith
,
R.
,
Y.
Chao
,
P.
Li
,
D.
Caron
,
B.
Jones
, and
G. S.
Sukhatme
,
2010
:
Planning and implementing trajectories for autonomous underwater vehicles to track evolving ocean processes based on predictions from a regional ocean model
.
Int. J. Robot. Res.
,
29
,
1475
1497
.
Sournia
,
A.
,
J. M.
Brylinski
,
S.
Dallot
,
P.
Le Corre
,
M.
Leveau
,
L.
Prieur
, and
C.
Froget
,
1990
:
Fronts hydrologiques au large des cotes francaises: Les sites-ateliers du programme Frontal
.
Oceanol. Acta
,
13
,
413
438
.
Stein
,
L. M.
,
1999
:
Interpolation of Spatial Data. Springer-Verlag, 246 pp
.
Stommel
,
H.
,
1989
:
The Slocum mission
.
Oceanography
,
2
,
22
25
.
Vignudelli
,
S.
,
P.
Cipollini
,
M.
Astraldi
,
G. P.
Gasparini
, and
G.
Manzella
,
2000
:
Integrated use of altimeter and in situ data for understanding the water exchanges between the Tyrrhenian and Ligurian Seas
.
J. Geophys. Res.
,
105
(
C8
),
19 649
19 663
.
Wallace
,
D.
,
2009
:
The surface ocean-lower atmosphere study (SOLAS). Global Change Newsletter, Vol. 73, Pacific Institute for Development, Environment, Security, 6–8
.
Zienkiewicz
,
O. C.
, and
R. L.
Taylor
,
1995
:
The finite element method. Basic Formulation and Linear Problems, O. C. Zienkiewicz and R. L. Taylor, Eds., McGraw-Hill, 648 pp
.