Automated feature tracking and vehicle navigation have the potential to facilitate autonomous surveys of ocean eddies by increasing sampling quality and/or decreasing operator workload. During an observational campaign in late 2013 and early 2014, methods for automated tracking were used to direct multiple ocean gliders during persistent surveys of a California Undercurrent eddy in Washington and British Columbia, Canada, coastal waters over a 3-month period. Glider observations of depth-averaged currents in the ocean’s upper kilometer and vertical separation of selected isopycnals were assimilated into a simple model describing eddy position, size, strength, and background flows using an extended Kalman filter. Though differing in detail from observations, results show the assumed eddy structure was sufficient to describe its essential characteristics and stably estimate eddy position through time. Forecast eddy positions and currents were used to select targets automatically to guide multiple gliders along transects through the eddy center as it translated. Transects performed under automated navigation had comparable or better straightness and distance from the eddy center when compared to navigation based on manual interpretation of the eddy scale and position. The tracking results were relatively insensitive to model choices at times when the eddy was well sampled, but they were more sensitive during sampling gaps and redundancies or rapid eddy translation. Overall, the results provide evidence that automated tracking and navigation are feasible with potential for widespread application in autonomous eddy surveys.
Mesoscale and submesoscale coherent eddies are ubiquitous in the upper ocean, account for a large fraction of oceanic kinetic energy at subinertial frequencies (Wunsch 2007; Ferrari and Wunsch 2009), and are responsible for significant redistribution of heat, freshwater, nutrients, and biota. A detailed understanding of the life cycle of eddies is critical for accurate prediction of where and when water and energy are transported (e.g., Danabasoglu et al. 1994), testing predictions of how eddies feed back onto the general circulation (e.g., Holland 1978), and for elucidating the role of eddies in the transfer of energy from large climatic input scales to microscales at which it is dissipated (e.g., Arbic et al. 2013), among many other topics. One approach that is used to investigate these processes is to persistently observe and monitor individual eddies.
Tracking of coherent eddies using gridded sea surface height (SSH) products from satellite-based radar altimetry has provided a dramatic increase in understanding of surface eddy frequency of occurrence, life-span, generation regions, trajectory characteristics, and general phenomenology (e.g., Chelton et al. 2007; Chaigneau et al. 2009; Chelton et al. 2011). However, the suite of gridded SSH fields presently in wide use has known limitations. An obvious drawback, the lack of resolution of subsurface fields, has been approached in an ensemble sense through analysis of Argo floats contained within eddies (Zhang et al. 2013; Lyman and Johnson 2015; Pegliasco et al. 2015). While providing extensive temporal and spatial coverage, gridded altimetry products attenuate the high-wavenumber portion of the mesoscale SSH spectrum (Chelton et al. 2011). This is problematic for the detection of subsurface-intensified eddies with weak surface expressions and/or horizontal extents that are of the same order as or smaller than the altimeter ground track spacing. Complementing gridded surface altimetry products through development of additional, widely applicable strategies for subsurface sampling of eddies is important to observational physical oceanography.
One method for complementing altimetry-based eddy life cycle studies is through buoyancy-driven autonomous underwater glider surveys. While gliders collect fewer types of observations than ships, and move more slowly (Rudnick and Cole 2011; Rudnick 2016), they have a much greater persistence than can be achieved from research vessel operations of comparable cost (Testor et al. 2010). During missions that can last 6 months or more (Pelland et al. 2016; Rudnick et al. 2016)—a duration equal to or greater than the lifetime of >90% of the eddies detected for 4 weeks or longer by Chelton et al. (2011)—gliders can obtain independent realizations of bulk eddy properties over time scales of days to weeks (Martin et al. 2009; Todd et al. 2009; Yu et al. 2017) while sampling at sufficient horizontal and vertical resolution to observe processes that moderate eddy evolution and decay (Pelland et al. 2013; Bosse et al. 2015, 2016, 2017; Pietri and Karstensen 2018).
Beginning October of 2013, an observational study was conducted off the Washington (United States) and British Columbia (Canada) coast, whose goal was to track and repeatedly survey a single California Undercurrent eddy (Cuddy; Huyer et al. 1998; Garfield et al. 1999; Collins et al. 2013; Pelland et al. 2013) using multiple Seaglider autonomous underwater gliders. Cuddies are one member of a large class of depth-intensified anticyclones found elsewhere with radii comparable to or smaller than the first baroclinic mode deformation radius (McWilliams 1985) and whose average surface signature and detectability in gridded altimetry products are subjects of active research (Stammer et al. 1991; Bashmachnikov and Carton 2012; Bashmachnikov et al. 2014). During these surveys, methods were developed to objectively track and forecast the subject eddy position, and to use these forecasts to automatically navigate the gliders. Objective eddy tracking and automatic navigation have the potential to facilitate autonomous surveys of eddies by increasing sampling quality or efficiency (Leonard et al. 2010; Flexas et al. 2018) while reducing operator workload.
This article describes the design and implementation of eddy surveys, provides an overview of the objective tracking and navigation methods and their results, and explores the sensitivity of the results to different assumptions or configurations. The results demonstrate the feasibility of automated tracking and navigation in autonomous eddy surveys with methods straightforwardly adaptable to other eddy surveys. Analyses of the eddy dynamics and evolution from these surveys are reported elsewhere (Steinberg et al. 2018, hereinafter SPE).
2. Field campaign
Cuddy surveys were performed using three Seaglider buoyancy-driven profiling vehicles (Eriksen et al. 2001) from the University of Washington School of Oceanography (designated SG189, SG194, and SG195, respectively; data from these deployments are available at https://accession.nodc.noaa.gov/0162344, https://accession.nodc.noaa.gov/0162349, and https://accession.nodc.noaa.gov/0162357, respectively). Additionally, data from a Seaglider (SG108) operated by the University of Washington Applied Physics Laboratory as part of the Northwest Association of Networked Ocean Observing Systems (NANOOS; http://www.nanoos.org) were used in aiding the search for a candidate eddy. Seagliders profiled from the surface to the shallower of either the bottom depth or 1000 m in a series of sawtooth dive–climb cycles; the School of Oceanography vehicles collected 1207 such cycles (2414 profiles) between initial deployment on 25 October 2013 and recovery of SG194 on 10 May 2014. The mean duration of all cycles with a maximum depth > 900 m was 8.65 h [1.40-h standard deviation (SD)] and the mean net displacement over ground was 6.31 km (1.98-km SD). The mean net horizontal speed through water—that is, net displacement through water estimated from a vehicle flight model (Eriksen et al. 2001), divided by cycle duration—was 21.40 cm s−1 (2.85 cm s−1 SD).
In this study, Seagliders collected samples of in situ conductivity, temperature T, pressure, optical backscatter and fluorescence, and dissolved oxygen using instrumentation and methods consistent with those described for other Seaglider deployments by Pelland et al. (2013, 2016). Seagliders also provide an estimate of depth-average current (DAC) over each dive–climb cycle, considered accurate to 1 cm s−1 (Eriksen et al. 2001; Pelland et al. 2016). Conductivity and temperature were measured using Sea-Bird Electronics (SBE) thermistor and conductivity cells mounted dorsally within a cylindrical duct, the flow through which is induced by vehicle motion rather than an electric pump (Perry et al. 2008; Pelland et al. 2013, 2016). Samples of conductivity, T, and pressure were used to estimate in situ salinity S, density ρ, potential temperature θ referenced to the sea surface, potential density referenced likewise, and spice (Flament 2002). The latter is used here as an approximate indicator of temperature or salinity anomaly along an isopycnal, with higher spice indicating warmer, saltier water for a given density. Vehicle sampling rates varied with depth, based on prescribed temporal sampling rates and vehicle vertical speeds. For School of Oceanography vehicles, the mean vertical spacing between samples of T and S was 0.98 m (0.50-m SD) shallower than 150-m depth, 2.00 m (0.81-m SD) between 150- and 300-m depth, and 3.83 m (1.47-m SD) deeper than 300 m. For SG108, the mean vertical spacing was 0.53 m (0.32-m SD) shallower than 50-m depth, 1.14 m (0.28 m SD) between 50- and 200-m depth, 2.34 m (0.33-m SD) between 200- and 500-m depth, and 3.28 m (0.86-m SD) deeper than 500 m.
b. Eddy acquisition
As part of the NANOOS network, SG108 performed repeat cross-shore transects between an inshore waypoint at 47°51.93′N, 125°10.27′W and an offshore waypoint at 47°N, 127°W during the period 24 September 2013–12 February 2014 (Fig. 1). During a transect performed between 15 and 24 October 2013, profiles from SG108 indicated two subsurface anticyclonic eddies carrying warm and saline water typical of Cuddies, one 20 km from the shelf break and another 107 km from the shelf break (Fig. 2). Two Seagliders (SG194 and SG195, respectively) were launched on 25 October 2013 at 47°N, 125°W and navigated toward the inshore candidate eddy, located at 47°51.00′N, 125°28.80′W. These vehicles traveled along the coast northwest of SG108’s transect line to approximately 48°15.00′N without detecting an eddy, but after being redirected to return southeastward, they encountered an eddy near the north rim of Quinault Canyon (Fig. 3).
On 20 November, SG194 first found evidence of increased isopycnal separation in the main pycnocline, elevated temperature and salinity along isopycnals, and cross-isobath offshore currents suggestive of anticyclonic flow (Figs. 3a,b). These features are all associated with typical Cuddies (Huyer et al. 1998; Collins et al. 2013; Pelland et al. 2013). SG195, which had proceeded farther to the southeast, was instructed to return north toward SG194, and it began sampling the target Cuddy on 24 November (Fig. 3c).
Surveys were conducted following Cuddy detection from 20 November until 25 February; further searches and possible surveys of a separate eddy, which are not discussed here, continued until the termination of the field experiment and transition to recovery operations on 21 April 2014. The survey strategy was to modify vehicle targets after each cycle in order to orient the two vehicles, SG194 and SG195, such that they performed repeat transits along orthogonal transects bisecting the eddy core and extending radially outward 20 km from its center. Manual tracking of the eddy was performed using an objective estimate of the eddy center, described below, as well as through subjective estimation of the eddy center position by examining spatial patterns of observed currents, isopycnal separation, and temperature/salinity anomaly along the core eddy isopycnal over several days in plan view. Near-real-time gridded sea surface height estimates were obtained from AVISO (http://www.aviso.altimetry.fr/) during the surveys, though the target eddy was not visible in these (data not shown). The eddy was tracked using only information retrieved by the four Seagliders involved in the field experiment (active survey vehicles shown in Fig. 1c).
The surveys during the period 20 November 2013–25 February 2014 were composed of three phases characterized by contrasting eddy translation speeds. The first phase, from 20 November to 6 January, consisted of poleward eddy movement along the continental slope until 3 December (markers 1 and 2 in Figs. 1a and 1b), followed by steady and slow movement offshore to the southwest (markers 2 and 3 in Figs. 1a and 1b). From 20 November 2013 to 6 January 2014, the eddy was estimated to have a net displacement of 65 km to the northwest over 47 days (Fig. 1a). The average daily movement during this period was 3.5 km.
Beginning 6 January, the results suggested that eddy translation had begun to accelerate in the poleward direction along the continental slope. SG194 was oriented to make sections in the along-slope direction of eddy propagation, while SG195 attempted to make cross-slope sections that were consistent in an eddy frame of reference (Fig. 1a). A third vehicle, SG189, was deployed to assist in the surveys on 16 January 2014. Automatic targeting, which will be described below, was initiated for SG195 during cycle 256 on 10 January 2014, for SG189 during cycle 55 on 30 January 2014, and for SG194 during cycle 309 on 1 February 2014 (Fig. 1c). Results suggested reduced eddy translation speed on 10 February 2014. During the period 6 January–10 February, the eddy was estimated to have traveled 255 km to the northwest over 35 days (markers 3 and 4 in Figs. 1a and 1b), with average daily movement of 8.0 km.
In the third phase, surveys were performed near the eddy’s estimated stopping position near the 2000-m isobath from 10 February onward (Fig. 1a). An electrical problem forced the withdrawal of SG195 on 21 February. On 25 February the eddy signal was no longer apparent in data returned from the remaining vehicles, and they were instructed to proceed offshore, which was judged to be the most likely direction of possible eddy escape based on the vehicle track history and the surrounding bathymetry. The estimated net eddy displacement from 10 to 25 February was 9 km to the west-northwest over 15 days (markers 4 and 5 in Figs. 1a and 1b), with average daily movement of 2.3 km. Further searches and surveys continued offshore until vehicles SG194 and SG189 were withdrawn because of persistent satellite communication difficulties on 3 and 21 April 2014, respectively.
3. Eddy model
An automated tracking method for the target Cuddy was developed with the goal of providing an up-to-date, objective estimate of the eddy center position and its translation velocity. The principle behind automated tracking was to form a simple kinematic model of the eddy that is described by a small number of parameters (e.g., eddy center position, propagation velocity, size, and strength), then estimate the time-varying values of these parameters using observations from the surveys. Two vortex kinematic models were tested. The first assumed a radial eddy velocity structure consistent with a Rankine vortex (Ide and Ghil 1998b; Kundu and Cohen 2008, section 3.11) and was used to forecast eddy position during surveys from 20 November to 24 December 2013. The kinematic model described below was found to improve upon the Rankine vortex assumption and was used thereafter.
The model assumed a geostrophic, baroclinic, axisymmetric vortex centered at time-varying Cartesian position . The potential density field was assumed to consist of a background profile that is a function of upward vertical coordinate z alone, and an eddy density anomaly that is a function of z and radial distance from the eddy center, :
and with a vertical structure function to be determined from the data. Here kg m−3 is a reference density, and the strength of the density anomaly is described by the dimensionless parameter .
The vertical structure was determined by computing the potential density anomaly within the eddy from an empirical estimate of the background profile . The background profile was first estimated on 3 January 2014 using Seaglider profiles collected from 18 December 2013 onward, when the vortex was well sampled during slow offshore propagation. Twenty Seaglider cycles sampling conditions consistent with ambient waters—that is, weak observed currents, high pycnocline stratification, and low values of spice and optical backscatter relative to those in the eddy core—were selected and an average over these was used as an estimate of .
A time series of potential density anomalies collected from SG194 as it traversed the eddy during the period 18 December 2013–3 January 2014 shows evidence of consistent anomalies of opposite sign surrounding the vortex core depth near 200 m (Fig. 4). The positive anomalies shallower than the core have a maximum strength of 0.2 kg m−3, while the deeper anomalies are 0.1 kg m−3 (Fig. 4). Vertical profiles of from both vehicles, normalized by the minimum value deeper than the vortex core (denoted as ), have a consistent shape, which supports the assumption of a separable spatial structure (Fig. 5a). There is an asymmetry in the vertical scale and amplitude of the density anomalies above and below the vortex core; plotting profiles in a stretched vertical coordinate , where is the background buoyancy frequency (Zhang et al. 2013, 2017), reduces the asymmetry in vertical scale (Fig. 5b). The median profile is approximated well by the first derivative of a function:
where α, β, γ, ω, and λ are constants; is the error function; and and correspond to the standard normal probability density and cumulative distribution functions (CDFs), respectively. The first term on the right-hand side of (3) is proportional to a “skew normal” probability density function (Azzalini 1985). The function (3) was fit to the integral , giving , , km, km, and . The derivative of this fit (red lines in Figs. 5a,b) was taken and, after substituting z coordinates for , was used as an estimate of . This form of does not have a dynamical basis; other possible approaches could include fitting a sum of low-order baroclinic modes (e.g., Cornuelle et al. 2000) or universal functional forms derived from Argo float data (Zhang et al. 2017). The profile , though not , was reestimated on 30 January 2014 as the vortex propagated close to the continental slope and the background stratification weakened.
If the vortex is assumed to be geostrophic, and that its swirl velocity is negligible at H = 1000-m depth, then the geostrophic azimuthal velocity , where , is given by
and 0–1000-m depth-average azimuthal velocity is given by
where g is gravitational acceleration, f is the Coriolis parameter, and is an amplitude term proportional to the strength parameter . By convention, A is negative for an anticyclonic eddy with , since ; through (5), this gives negative (clockwise) azimuthal depth-averaged currents. In all subsequent notation, the angle brackets around velocity vectors or components are dropped, since all velocities referred to in this article are 0–1000-m depth averages.
The aforementioned system includes the parameters R, A, , and , which collectively describe the structure of eddy-induced currents and density anomalies. Idealized models of Mediterranean outflow eddies, which are configured similarly to Cuddies, suggest that in the presence of energetic boundary currents such as near the Washington continental slope, eddy translation caused by self-propagation is likely to be overwhelmed by advection within background currents (Dewar and Meng 1995). This is consistent with strong along-slope translation velocities inferred by repeat observations of eddies in Pelland et al. (2013). Therefore, it was further assumed that eddy translation was due to advection by a background horizontally uniform barotropic current , such that
and that observed depth-average currents were a superposition of eddy currents and the background values
where are azimuthal currents resulting from the eddy as estimated from (5), converted to Cartesian coordinates.
4. Estimation, prediction, and navigation
The translating eddy system is described by the parameters R, A, , , U, and V, which together compose the elements of a state vector , where a superscript T indicates a transpose. Seaglider observations were used to estimate using the discrete-time extended Kalman filter (EKF), a widely used tool in estimation and control of dynamical systems, including in atmospheric or oceanographic applications for data assimilation in general circulation models (e.g., Evensen 1992; Gauthier et al. 1993; Fukumori et al. 1993; Wunsch 1996; Fukumori et al. 1999). The EKF has also been tested as a method for tracking idealized systems of vortices by Ide and Ghil (1998a,b). The EKF was used in this application because it is a method for sequential estimation of the state of a dynamical system, when either the time derivative of the system state or the observations of the system are a nonlinear function of the system state. Here, while is a linear function of —all parameters are modeled as persisting in time except the center location, which is related to translation velocity by (6)—observations of both velocity and density anomaly from Seagliders are nonlinear functions of the state parameters because of the Gaussian radial structure assumed in (2) and (5).
It is assumed that in the time interval between two sets of observations (available at times and ), the state vector evolves according to
where is the transition matrix; and is the “process noise” vector, which represents the residual evolution of eddy state parameters as a result of unresolved processes not explicitly included in the eddy model. The process noise vector is assumed to be composed of zero-mean, normally distributed, time-uncorrelated random perturbations (white noise) with covariance matrix , where the operator indicates the expectation of a random variable (Gelb 1974, p. 30). The specification of the elements of is described below [section 4c; (18)]. The transition matrix expresses persistence of the eddy size and strength together with advection by the background flow.
Sets of observations of the system (currents, isopycnal separation; section 4b) at each time are arranged in the vector and are assumed to obey
where is an vector of observations that would be expected in absence of any observation noise. This vector is a nonlinear function of system state . The number of observations may vary at each time. The observation noise—residuals resulting from sensor error and unresolved processes—is represented by the vector . It is also assumed to be normally distributed white noise, with covariance matrix independent of the process noise. For element i of , which corresponds to an observation collected at position (, ) and time , the corresponding element of is the expected component of current [from (5) and (7)] or isopycnal separation [from (1) and (2)] based on the eddy state at . As an example, if observation i corresponds to the zonal component of a DAC observation , then element i of is , where
is the radius of observation i from the eddy at ; is the x component of a unit vector in the direction of positive azimuthal currents; while , , and are the eddy strength, eddy characteristic radius, and zonal component of background currents, respectively, as above.
In this study, for each forward run of the EKF, a starting guess was provided (Table 1), along with an estimate of the error covariance of this starting guess, , described below. Subsequent to the time of the starting guess, Seagliders would surface to transmit profile and engineering data after each dive–climb cycle. Following a successful transmission, raw data were processed (section 2) to yield estimates of T, S, ρ, and , available within a few minutes of cycle completion. A new estimate of the state vector was composed each time a successful profile data transmission occurred. The interval between transmissions was on average 4 h but was as long as 10 h and as short as a few minutes during the surveys.
The vector of differences between and the vector , which represents the expected observations, is denoted as the innovation :
Here and is the previous best estimate made after the observation window . Note from (12) that since the expected observations take into account the eddy propagation with time, and expresses linear translation and persistence of all other properties, for practical purposes, .
An updated estimate of the state vector, , was then produced by adjusting the estimate based on the product of the innovation and the gain matrix :
where the gain is given by (Gelb 1974, 186–187)
with superscript −1 indicating a matrix inverse. In (15), is an matrix of partial derivatives of the expected observations with respect to the forecasted state parameters, here evaluated using centered finite differences of the forecasted state vector elements. The matrix
is the state error covariance estimate, forecasted at from its previous updated estimate . An updated estimate of the state error covariance at , , was then computed according to (Gauthier et al. 1993)
where is the 6 × 6 identity matrix.
In the limit where the nonlinearities in the observation operator in (9) vanish, the gain [(15)] minimizes the expected square error of the estimated state, under the assumptions made in (8) and (9), and the additional assumption that and errors in are uncorrelated (Gelb 1974; Gauthier et al. 1993; Brown and Hwang 1997). Where the nonlinearities are nonzero, the EKF instead provides a first-order approximation to the optimal state estimate update (Anderson and Moore 1979; Ide and Ghil 1998a).
The observations used in were the zonal and meridional components of DAC on each Seaglider cycle, used throughout the surveys, and the vertical separation between two chosen isopycnals on both the descent and ascent profiles of each cycle, used following 7 January 2014. The DAC components measure eddy swirl velocity and background currents, while the isopycnal separation is a measure of eddy density anomaly. Spice along a core isopycnal was also considered but was rejected, since the background conditions were difficult to define because of ambient mesoscale noise and strong cross-shore and alongshore structure (Pelland et al. 2013). The horizontal gradient of spice is a potentially useful observation type (Fig. 2) that could be considered in future studies of eddies with strong water mass contrasts from their surroundings.
The isopycnal vertical separation was computed on each profile between the σθ (≡ ρθ − 1000 kg m−3) = 26.5 kg m−3 and σθ = 26.8 kg m−3 isopycnals () . These were selected based on the expectation that they would have the greatest upward and downward vertical displacement, respectively, from their background levels given the estimated profile and empirical eddy density anomaly vertical structure function . These isopycnals were not altered following reestimation of on 30 January 2014. For EKF state estimates made after 30 January, values collected after 15 January were compared to those expected based on the second estimate of .
For each eddy state estimate, observational data from the most recent six cycles of each vehicle were included in , provided that these cycles were not performed more than 56 h prior. That is, the subscript k on this vector refers not to observations collected precisely at , but rather a set of observations collected at different locations and times up to 56 h prior to . This 56-h time window was chosen because it was assumed that using many more observations than state vector elements would more heavily constrain the state update at each step, which otherwise depends on the poorly known measurement and process noise variances. The trade-off is that, since each vehicle performed many dive–climb cycles within a 56-h period and the eddy state was updated after each cycle, observations from any single cycle were assimilated in multiple windows. As a result, each windowed set of observations was not independent. This choice is discussed further in section 6.
c. Model noise
Beyond model stability, an additional significant challenge in the implementation of the EKF is that the character of the model and process noise has a large influence on the performance of the filter, though it is often the case that noise values are not well known a priori (Gauthier et al. 1993). If the measurement noise dominates over the process noise , then the denominator becomes large relative to the numerator in (15) and the Kalman gain is small, hence the estimated state is not very sensitive to new observations. Conversely, the filter can become oversensitive to small errors in the observations if the measurement noise is assumed to be small relative to the process noise.
In this application, we assumed that the matrices and were diagonal, with elements along the main diagonal equal to a noise variance specific to the type of state vector element or observation for that row. For the process noise, it was assumed that
where km, m s−1, m, and m2 s−1. For the measurement noise covariance matrix , whose size changed at each time step according to the number of observations within each window, diagonal elements were m s−1)2 for rows corresponding to DAC components and m)2 for rows corresponding to isopycnal separation. Off-diagonal elements of were set to zero.
It should be emphasized that these noise variances were chosen experimentally by selecting values that allowed a compromise between forecast stability during relatively steady eddy translation and responsiveness to rapid accelerations; information about the statistical properties of continuous Cuddy translation and evolution was not available. In the selection process, the measurement noise standard deviations were restricted to be greater than the absolute sampling accuracy of Seagliders (section 2a), because in the EKF, measurement noise must take into account both error caused by instrument accuracy and variance in the observations caused by processes other than the eddy-background system (e.g., tides, inertial oscillations, internal waves), which for the purposes of the model represent measurement noise. The assumed noise variances described here were used in eddy forecasts from 15 January 2014 onward.
Seagliders in this and other missions are directed using “targets” files (formatted text files containing a list of navigational targets through which the vehicle is instructed to proceed) that the vehicle downloads during its surface communication sessions with the operator base station computer via Iridium satellite telemetry. For manual navigation in this study, it was often the case that operators were required to update these files multiple times each day for each vehicle. To automate navigation an algorithm was devised using the EKF eddy state output to reduce operator workload and to objectively choose targets. Automated navigation was used for 58 (124) cycles of SG189, 59 (20) cycles of SG194, and 121 (0) cycles of SG195 before (after) 25 February (Fig. 1).
The goal of the algorithm was to compute a target that would guide the vehicle along a transect across the eddy. The first step of the algorithm was to forecast the next surfacing location of the vehicle based on the latest estimate of eddy position, size, strength, and background currents. The target was then chosen based on one of two cases (Fig. 6):
Continuing an existing transect of the eddy (Fig. 6a). Transects were intended to follow a uniform course in an eddy-centric polar coordinate system. Each target was chosen such that in absence of error, estimated currents and vehicle displacement through water would result in future surfacing along the desired transect.
Initiating a new transect of the eddy (Fig. 6b). This occurred if the vehicle was estimated to surface beyond an operator-specified turnaround radius from the eddy center (). The azimuth of this new transect was defined based on the bearing from the forecasted surfacing location to the eddy center. The value of was originally set at 15 km and modified to 10 km on 13 January.
Following case 1, as the vehicle proceeded toward the eddy center, the algorithm would allow the transect azimuth to be redefined based on each new surfacing, until the vehicle was forecasted to surface within the characteristic R, after which it was held fixed. If the forecasted eddy currents were too strong for the vehicle to surface on the desired transect, a scenario not shown in Fig. 6, then the target was set such that the vehicle would choose a heading in the desired transect direction in an attempt to stem the currents.
The EKF and navigation algorithms were implemented on the University of Washington School of Oceanography Seaglider base station desktop computer using the Octave programming language (https://www.gnu.org/software/octave/) within a function that ran automatically following each upload and successful processing of Seaglider data,1 and which could also be run on command from an operator at any time. The code package used during the field experiment is available upon request from the corresponding author.
e. Configurations and reference run
As noted above, the model configuration used in eddy forecasts changed over time. The EKF was first started on 20 November 2013 under the assumption of a Rankine vortex radial structure, and only observations of were assimilated into the model. On 13 December 2013, the eddy state vector time history was reestimated based on a new starting guess from 8 December 2013. This was done because the eddy center position estimated by the model was appearing to diverge from the subjectively estimated position at that time. On 24 December, the state vector was reestimated from 8 December but under the assumption of a Gaussian vortex structure and a smaller starting estimate of R. On 6 January, the state vector was again reestimated from 8 December but with isopycnal separation as an added observation type. The starting guesses for these various model configurations are described in Table 1. In all cases, the assumed starting error covariance was set equal to a 6 × 6 zero matrix; in practice, the EKF generally converges on an equilibrium in the error covariance matrix within a few time steps regardless of the starting guess. Configuration changes were prepared offline before rerunning the EKF function while all vehicles were submerged, hence configuration changes did not result in interruptions to tracking.
After completion of the observational campaign, a post hoc forward run of the EKF was performed from 20 November 2013 to 25 February 2014 using the final eddy model, observations, process noise, and measurement noise described above. This run, denoted the “reference” forward run, differs from what was used in the original surveys, in which the final forward run was started with an initial guess on 8 December 2013 (Table 1).
Note that for an estimate at time , this reference run provides an estimate of the eddy system state that reflects observations collected prior to , but not after. In contrast to a forward run (“filtering”), forming an estimate that takes into account data both before and after is denoted the “smoothing” problem. Smoothed estimates of the eddy state vector were also computed in this study, using the reference forward run as a basis, to provide an estimate of the eddy position versus time that includes all available data. The smoother equations and implementation are described in the appendix.
a. Eddy track, structure, and sampling
In what follows, references to the “forecasts” indicate the best estimates of the eddy properties that were available during the observational campaign, using the model configurations summarized in the first four columns of Table 1. The “smoother” refers to the smoother estimate of the eddy track made using the final version of the model, as described above (section 4e) and in the rightmost column in Table 1.
The forecast and smoother eddy center tracks are shown for comparison in plan view in Fig. 7a, while the differences in zonal and meridional position versus time between the two are shown in Figs. 7b and 7c, respectively. Following 24 December 2013, when the forecasts were modified to assume a Gaussian rather than Rankine radial velocity structure, the two tracks are generally within a few kilometers of one another and the difference in position is within the estimated 95% confidence bounds on the error of the forecast track. The median position difference between the two is 1.7 km. The greatest deviations occur on 11 December 2013 (13.2 km) and 17 January 2014 (12.7 km). These periods will be discussed in more detail below; the results suggest that the vehicle sampling characteristics and, in the latter case, rapid eddy translation (Fig. 7d) resulted in increased uncertainty in the eddy position.
The estimates of the background flow are in qualitative agreement between the forecasts and smoother, though differ in detail during some periods, for example, 20 November–1 December and 8–15 December 2013 (Fig. 7d). For both of these periods, the forecast track assumed a Rankine vortex structure. Estimated 95% confidence bounds on the background flow are of order ±1 cm s−1 in either the forecasts or the smoother.
Vehicle zonal and meridional position relative to the forecast track over 20 November 2013–25 February 2014 is shown in Fig. 8a. A total of 1266 profiles were collected within 30 km of the forecast eddy center during this period, 559 (44.2%) of which were within 10 km. Here, the location of each slantwise profile is defined as the radius at which the profile crossed the eddy core depth. Results are similar when comparing vehicle position to the smoother track (Fig. 8b). Vehicles collected 1280 profiles within 30 km of the smoother position, of which 560 (43.8%) were within 10 km.
Observed 0–1000-m DAC relative to the smoother eddy center position is shown in Fig. 8c. A clear signal of anticyclonic azimuthal flow around the eddy center is evident, though with an apparent asymmetry in which flow is weaker in the northeast quadrant. This asymmetry likely reflects a superposition of eddy currents and the net northwestward background flow associated with poleward eddy translation along the continental slope (Fig. 7d). Removing the smoother-estimated background flow from each observation (Fig. 8d) yields a clearer eddy signal.
Figure 8e shows normalized observations of azimuthal currents versus normalized radius . Here is the azimuthal component of the currents in Fig. 8d, divided by the smoother-estimated peak currents at that time, [cf. (5)]. Consistent with Fig. 8d, though there is substantial noise in , anticyclonic flow of the correct magnitude and approximate radial structure is evident. If the eddy radial structure were exactly as assumed in (5) and all estimates were perfect, then the normalization chosen in this figure would collapse all observations along the red curve in Fig. 8e. An average radial profile of in bins of (black curve in Fig. 8e, dashed curves ±1 SD) is in good agreement with the assumed radial structure to at least the radius of maximum velocity. Beyond this radius, the composite exhibits a weaker radial decay than assumed.
Figure 8f shows normalized isopycnal separation versus . Observed is in this case obtained by subtracting the background and dividing by the maximum increase in from the background that would be expected at the eddy center. The red curve illustrates a Gaussian radial decay , which is the approximate radial structure of nondimensionalized isopycnal separation that would be expected2 from the eddy model [(2)]. The composite profile of has a similar radial decay, though it is weaker than expected for .
The estimated maximum strength of depth-average eddy currents was on average 6.33 cm s−1 in the forecasts, while the radius of maximum eddy currents was on average 10.01 km. Composing averages of these metrics, as opposed to the quantities A or R, allows averaging across the entire forecast period, during which both Rankine and Gaussian vortex parameterizations were used (Table 1). By comparison, the average values from the Gaussian smoother are 5.55 cm s−1 and 7.99 km, respectively. Since the focus of this study is on the tracking and forecasting of eddy position, these properties are not discussed further here. For a more detailed analysis of eddy kinematic and dynamical properties, and their evolution with time, see SPE.
b. Automated navigation
The vehicle survey tracks relative to the smoother eddy center obtained under automated navigation are qualitatively similar to those obtained via manual navigation for the period prior to 25 February (Fig. 9a). During this period, 68% of profiles performed under automated navigation were within 10 km of the smoother eddy center. The quality of the eddy surveys is compared between manual and automatic navigation here using two criteria: the radial distribution of profiles and the character of individual transects of the eddy. Overall, the radial distribution of profiles collected under automated navigation was offset toward the eddy center, with a significantly smaller mean radius and interquartile range (IQR; Fig. 9b; Table 2) than the distribution under manual navigation.
To identify individual transects of the eddy for each vehicle, the vehicle radial distance from the smoother eddy center versus time was computed. Transects were defined as periods between successive maxima in this distance, between which the vehicle closed to within at least 10 km of the eddy core; 98 total transects were identified, with 51 of these performed entirely under automated navigation. Transects are here evaluated using two metrics: one related to deviation from a straight line (“deviation”) and one related to how closely each transect comes to bisecting the eddy (“bias”).
Transect bias and deviation were evaluated using profile locations in an eddy-centric coordinate system (see example in Fig. 10a). The main axis of each transect was defined as a line that minimizes the sum of square residuals along a direction orthogonal to that line; bias was defined as the minimum distance between this line and the eddy center, while deviation was evaluated as the root-mean-square (RMS) residuals. Overall, average bias and deviation were significantly less under automated navigation than under manual navigation (Figs. 10b,c; Table 2).
Automated transects were significantly shorter than manual (Table 2). Both types were purposefully shortened during periods of rapid eddy translation (Fig. 10d). There is only a limited period in which both 14 automated and 10 manual transects were collected simultaneously (10–30 January 2014; Fig. 10d). During this period, automated transects had a significantly lower bias but a similar deviation (Table 2). Transects were of comparable length, and the radial distribution of profiles were also similar, under either navigation type during the overlap period (Table 2).
Using the data obtained from the eddy surveys, it is possible to evaluate the sensitivity of the eddy forecasts to alterations in the assumed model or observational data. Since the focus of this study is tracking the eddy position, sensitivity is evaluated by altering the model or observations from those used in the reference forward run, and examining the differences in position from this run over time. Note that in a real survey, changes in eddy position estimates may lead to differences in sampling, an effect that is not considered here.
The measurement and process noise levels were not constrained by observational data and were selected experimentally (section 4c). To evaluate sensitivity to the assumed noise levels, for each of the six unique variance parameters (measurement: , ; process: , , , ) four additional post hoc forward runs were performed, in which the chosen parameter was increased and decreased by factors of 2 and 10. The minimum/maximum, median, and first/third quartiles of the eddy position differences from the reference forward run—denoted as —were then computed. The distribution of the values for each run are shown in box-and-whisker plots in Fig. 11a. All runs in which the noise variances were altered by a factor of 2 had median from the reference run of less than 1 km (Fig. 11a). Median was greater for runs in which variances were altered by a factor of 10 and up to 5.2 km in the run in which was increased by 10. In plan view, the trajectory of this run noticeably diverges from the reference run, in contrast to the remaining runs (Fig. 11b). Forward runs were in general more sensitive to modification of the measurement noise values than the process noise, with the exception of , the noise variance associated with perturbations to the background velocity (Fig. 11a).
The values of were not uniform versus time in each run and exhibited increased run-to-run spread at certain times (Figs. 11b,c). The median value versus time of across the 24 runs, a measure of run-to-run spread, was positively correlated with the 95% bound on the size of eddy position uncertainty in the reference forward run (). Two periods in particular are associated with increased run-to-run spread/uncertainty: 9–11 December 2013 and 16–19 January 2014 (Fig. 11c). The divergence of the run occurs following the latter period (Fig. 11c). In the former period, the forecasts under the Rankine model began to diverge from the subjectively estimated eddy center and were restarted on 13 December. The of the last Rankine run prior to the restart is shown in blue in Fig. 11c.
Additional forward runs were performed in which the observation types were limited to either or (Fig. 12a), or the state vector was reduced to four (position, background velocity) or two elements (position only, assuming zero background flow; Fig. 12b). In the latter two cases, the eddy strength A and R parameters were assumed constant and equal to the starting values from the reference forward run (Table 1). Differences in position from the reference run were overall more sensitive to these observation or state vector changes than the noise variances (Figs. 12a,b). Runs that assimilated only (Fig. 12a) or omitted background flow from the model (Fig. 12b) noticeably diverged from the reference run during the period of strong poleward propagation, likely because of the inability of the model to account for background flow in the observations (latter case) or the inability of the model to distinguish changes in background flow from eddy strength/position changes, because of a lack of isopycnal separation information (former case).
The agreement in forecast and smoother tracks, and the appearance of robust anticyclonic depth-averaged flow in an eddy-centric coordinate system, is evidence that the forecasts yielded realistic estimates of the eddy center position and that the surveys effectively captured its bulk structure. In contrast to other glider eddy surveys (e.g., Martin et al. 2009, their Fig. 7), the vehicle tracks in this study show modest, if any, distortion by the eddy currents. Here, the average vehicle horizontal speeds through water (section 2a) were almost 4 times as large as the average maximum 0–1000-m eddy currents estimated by the smoother (section 5a).
Rather than strong eddy currents, the main challenge to accurate survey of the eddy was its small size and rapid acceleration and translation, which motivated keeping vehicles near the eddy interior to obtain information about the center location from DAC and isopycnal separation. As a result, this inner region was densely sampled, with the trade-off that temporal variability in the background conditions was sparsely resolved. This reflects one of the principal limitations of buoyancy-driven gliders in eddy surveys, which is vehicle horizontal speed and the consequent multiday period required to obtain a single transect across the eddy. Speed and profiling rate were conservatively chosen in this application in order to maximize vehicle endurance, and could be modified in future studies. Increasing vehicle speed could reduce the time between transects of the eddy; its impact on tracking accuracy and stability using the EKF or other algorithms warrants further investigation. Sampling could conceivably be modified as eddy or background characteristics changed. The rapid translation of the target eddy in this study was not expected based on previously published Cuddy translation speeds (Collins et al. 2013, their Table 1). Estimates of ambient currents from SG108 support the hypothesis that the eddy movement was due to background flows (SPE).
Although differing in certain details in comparison to the observations, the very simple kinematic model chosen in this study was sufficient to capture the essential characteristics of the system and stably estimate the eddy position. One simplification was the choice of a geostrophic, rather than gradient wind, horizontal momentum balance, which for a given pressure gradient leads to an underestimate of the currents in anticylonic eddies (Elliott and Sanford 1986). While more realistic, a gradient-wind balance imposes narrower limits on the physically allowable combinations of eddy size and strength (McWilliams 1985); from a practical standpoint, with no other modifications, the filter could choose combinations of A and R that produce complex-valued eddy currents because of the quadratic gradient-wind velocity solutions (McWilliams 1985; Pelland et al. 2013). Rather than attempt to accommodate these restrictions by further modifying the filter, we elected to proceed under the geostrophic assumption, which was validated by the results. This assumption may, however, contribute to the observed misfit in isopycnal separation, in addition to the coarse time estimation of and background stratification during eddy tracking. In the analyses of SPE, the assumption of slowly varying vertical structure and background stratification is relaxed.
The eddy signal was eventually lost over the continental slope, where the assumptions of the kinematic model, especially those of background spatial homogeneity in the density field and currents, are likely to be least valid. The possibility that the eddy was ultimately destroyed because of background shear or interaction with topography also cannot be ruled out (Brickman and Ruddick 1990; Richardson and Tychensky 1998; Cenedese 2002; Torres and Gomez-Valdes 2017).
Automated transects were initiated during a period of rapid eddy translation and were directed to be short, in order to decrease the interval between transects and lower the likelihood of losing the target eddy. This resulted in a more confined radial distribution of sampling under automated navigation. It is possible that the overall shorter length of automated transects contributed to their reduced bias and deviation in comparison to manual transects, though transect length was not a predictor of bias and only a weak predictor of deviation. Within the very short overlap period, the automated transects performed at least as well as the manual transects by these metrics, and had a similar radial distribution of profiles. Straightness and bias were used as simple evaluation criteria because these were the goals of manual navigation in the field experiment described here. A more comprehensive evaluation of automated navigation could include its impact on tracking accuracy and stability, the estimation of eddy properties, and overall piloting workload. These effects are likely application specific, and their evaluation would require extensive trials or simulation of automated versus manual navigation schemes using observational or synthetic data.
Two distinct periods were apparent in which the eddy position was more uncertain than at other points during tracking (section 5c). During the first peak in uncertainty (Fig. 13a), both survey vehicles had been outside the eddy interior and were opposite the direction of propagation for the preceding 2 days. In the second peak (Fig. 13b), observations from both vehicles were concentrated in a comparatively small region during a period of rapid eddy translation. In contrast, during periods of relatively low uncertainty (Figs. 13c,d), the vehicles obtained samples over a wide region of both the eddy interior and exterior, along nearly orthogonal lines, providing a stronger constraint on the eddy center position, even during rapid eddy translation in early February (Fig. 13d). Though a complete investigation of the error and stability characteristics of the nonlinear eddy estimation algorithm is beyond the scope of this study, these examples suggest that crossing the eddy in multiple directions and maintaining continuous sampling of the interior region by at least one vehicle may reduce tracking uncertainty.
As discussed above (section 4b), an important stipulation made in eddy tracking was the use of overlapping 56-h windows when assimilating observations into the eddy state estimate, which violates the implicit assumption that the observations at are of the eddy process between and (an interval generally much less than 56 h; section 4a). An additional forward run was performed in which only a single cycle was assimilated at a time while holding all other properties equal to the reference run. This ensures that each set of observations is independent and represents new information since the last eddy state estimate. In this “single cycle” run, the observation vector is , where and indicate the isopycnal separation on the descent and ascent profiles, respectively, for cycle k. This run has a qualitatively similar trajectory to the reference forward run (Fig. 14), with a median difference of 3.2 km but a maximum difference of 19.4 km.
Aside from overlapping observation windows, another key stipulation in the eddy model and assimilation scheme was to assume that process noise variance was constant between updates, despite unequal time between these updates, and that there was no correlation between different process noise elements. A more accurate method of forecasting the state error covariance rather than (16) would have been to integrate the continuous-time error propagation equation between updates, which for the linear state process considered here is
where and is a continuous-time white noise process such that , where Δ(τ) is the Dirac delta function (Gelb 1974, chapter 3). This approach is sometimes known as a “continuous discrete” filter, in which a continuous process is observed at discrete intervals (Gelb 1974). A continuous-discrete run was performed while assimilating only a single cycle, and with , where h is the average time between update steps in the reference forward run and is as described in (18). This run is very similar to the single-cycle discrete run (Fig. 14), with a median of 3.1 km and a maximum of 21.4 km. These runs demonstrate that it would have been possible to generate a qualitatively similar eddy trajectory under automated tracking without the assimilation of multiple cycles, or the discrete-time approximation. However, since these analyses were performed after the completion of the surveys, their potential impact on the success of eddy tracking is not known.
Filter “consistency” refers to the degree to which the assumed measurement and process model and noise values are consistent in a statistical sense with the observations. The EKF can be quantitatively tested for consistency using the same tests that are applied to the linear Kalman filter (Bar-Shalom et al. 2001, section 10.3), which use the measurement innovations . Eddy state forecasts or the reference forward run were not evaluated for consistency in this study because the tests assume that the observation vector is of a consistent dimension throughout the filter run, and in the forecasts and reference run it is not (section 4b). The additional single-cycle runs, where throughout the run, were tested and were found to fail the statistical tests for consistency described by Bar-Shalom et al. (2001, section 5.4); these results are omitted here for brevity. As the survey results show, an EKF can provide useful information even in cases where some of its underlying assumptions are likely violated. Further exploration of the modifications to the filter necessary for consistency is beyond the scope of this study, though we note that building a consistent filter in future applications would benefit from a more thorough a priori characterization of the statistical properties of instantaneous eddy movement, which are presently poorly known. Where feasible, these could be addressed by altimetry observations of other eddies in the region of interest. During surveys, noise variances could also be estimated in parallel with the eddy state vector using an adaptive estimation algorithm (Brown and Hwang 1997).
In this study a California Undercurrent eddy was tracked for several weeks using only near-real-time information retrieved from four Seagliders. The success of the EKF forecasts and smoother in describing the eddy trajectory based on a simple kinematic model supports the idea that the use of such an objective method is a tractable approach for facilitating persistent in situ eddy surveys from gliders. A simple automated navigation routine based on this model successfully guided three vehicles for portions of the eddy surveys, leading to transects across the eddy center that were of comparable or better precision than those resulting from manual navigation, as judged by straightness and distance from the eddy center. Although this study does not provide a comprehensive test of automated navigation, the results provide evidence that it is a feasible strategy that can be implemented while at least maintaining overall sampling quality.
For the model used in this study, the choice of the unknown noise variances did not significantly influence the position estimates when the eddy was well sampled, but it did have a much larger effect during sampling gaps or redundancies. In a forecasting situation, it is possible this could have important implications for model stability and tracking success. Eddy position estimates in this study were sensitive to the removal of isopycnal separation observations, or to the removal of background flow from the eddy model, which in this case was a large component of the signal. In buoyancy-driven glider surveys, vehicle speed limits the rate at which independent realizations of the eddy are obtained; the use of multiple observation types helps to offset this limitation. After surveys were completed, modifications to the data assimilation scheme were made to correct known deficiencies in the original forecasts, which resulted in moderate changes to the estimated eddy position through time.
The methods used here have the potential for significant refinement or expansion. As examples, higher-order measurement or state evolution terms, or navigation schemes informed by objective sampling criteria (Eichhorn 2015; Zamuda et al. 2016; Lermusiaux et al. 2017), could be incorporated. The EKF could be adapted to assimilate data from additional sources, such as ships, profiling floats, or drifters, and in principle could also be paired with a more complete dynamical model of the eddy system (Flexas et al. 2018). Such a model could be used for the integrated purposes of both tracking and later analysis, which is in contrast to the approach here of choosing a tracking model that was as simple as possible, then refining the analysis after the completion of the surveys (SPE).
We emphasize that this study serves as a proof of concept rather than as a general exploration of the optimal strategy for, or limits of, persistent surveys of eddies using gliders. The methods described here were developed to track a small fast-moving vortex in a dynamic near-coastal environment using sampling optimized for endurance; it is likely that other surveys may face very different constraints or challenges. The choices of model complexity, noise variances, observation type, the number of vehicles, and their sampling strategy may all have important impacts on tracking success, which will also likely depend on the eddy and region of interest. Performing survey simulations using output from realistic regional numerical models (e.g., Kurapov et al. 2017; Flexas et al. 2018) is one possible approach to further evaluation of methods for tracking and automated navigation, which would allow repeated experiments over a wide range of eddy and background conditions.
The National Science Foundation’s Division of Ocean Sciences supported this work under Grant OCE-1153980. This work was completed while NAP was a postdoctoral investigator at the NOAA/Alaska Fisheries Science Center/Marine Mammal Laboratory and Joint Institute for the Study of the Atmosphere and Ocean. This study would have been considerably more difficult without the efforts of Craig Lee, Jason Gobat, Ben Jokinen, and Geoff Shilling in providing and operating Seagliders for use in NANOOS. We are grateful for their work and for sharing data from SG108. The authors are deeply indebted to Kirk O’Donnell and Bill Fredericks for their resourcefulness, dedication, and endurance of difficult conditions, which made the field operations in this experiment possible. Thanks to Tom Young and the crew of the F/V Tommycod for their professionalism and capability. We are also grateful to Bill Crawford, Doug Yelland, Marie Robert, and the officers and crew of the CCGS John P. Tully and U.S. Coast Guard Station Grays Harbor for their generous and timely assistance in recovering Seagliders used in this field experiment. Sarah Webster provided helpful discussion regarding Kalman filtering. We thank three anonymous reviewers and editor Dr. Denis Volkov for their thoughtful comments and editorial judgment, all of which helped to significantly clarify and strengthen a previous version of this manuscript.
These equations are a backward recursion in which the smoother estimate at depends on the forward-run estimate at , ; its error covariance matrix ; and the smoother estimate at , ; along with its error covariance matrix . When applying the smoother, the recursion is started at the second-to-last estimation time , where K is the total number of vehicle cycles used in eddy estimation, and with , since there are no observations after the final state estimate and therefore the smoother and filter estimates are equal there. To produce the smoother estimate, (A1)–(A3) were applied successively backward starting with the second-to-last state estimate from the reference forward run.
Current affiliation: NOAA/Alaska Fisheries Science Center, Seattle, Washington.
This article has a companion article which can be found at http://journals.ametsoc.org/doi/abs/10.1175/JPO-D-18-0033.1
When a Seaglider surfaces, it downloads any available targets file prior to uploading the data from its most recent dive–climb cycle to the base station. Therefore, if the vehicle is not directed to hold on the surface, it will submerge and initiate its next cycle before the EKF and autotargeting are updated based on its most recent data. This introduced a latency in autotargeting in this study: data returned from a vehicle on cycle j were not used in its navigation until cycle . This compromise was considered acceptable because in Seaglider piloting, it is not generally advisable to issue automatic commands for the vehicle to hold on the surface and then resume diving, as the latter may override important recovery procedures that activate when a vehicle detects a system fault.
Because the depth of a given isopycnal is a nonlinear function of eddy density anomaly, the expected radial structure of isopycnal separation differs slightly from this form at each observation point, based on that point’s combination of background stratification, eddy strength, and eddy radius.
Here for a two-sided Student’s t test of the correlation coefficient at 95% confidence with degrees of freedom (Panofsky and Brier 1968, section 4.6).
Here and .