## Abstract

Idealized simulations of autonomous underwater glider sampling along sawtooth vertical–horizontal paths are carried out in two high-resolution ocean numerical models to explore the accuracy of isopycnal vertical displacement and geostrophic velocity profile estimates. The effects of glider flight speed, sampling pattern geometry, and measurement noise on velocity profile accuracy are explored to interpret recent full-ocean-depth Deepglider observations and provide sampling recommendations for glider missions. The average magnitude of velocity error profiles, defined as the difference between simulated glider-sampled geostrophic velocity profile estimates and model velocity profiles averaged over the spatial and temporal extent of corresponding simulated glider paths, is less than 0.02 m s^{−1} over most of the water column. This accuracy and the accuracy of glider geostrophic shear profile estimates are dependent on the ratio of mesoscale eddy to internal wave velocity amplitude. Projection of normal modes onto full-depth vertical profiles of model and simulated glider isopycnal vertical displacement and geostrophic velocity demonstrates that gliders are capable of resolving barotropic and baroclinic structure through at least the eighth baroclinic mode.

## 1. Introduction

Since their development nearly 20 years ago, buoyancy-driven autonomous underwater vehicles (gliders) have been increasingly used to collect sustained observations of ocean water property and velocity fields. Their successful use in diverse environmental conditions achieving various mission objectives is the result of increased accessibility to the platform, as well as improved endurance and reliability (Rudnick et al. 2016). While gliders provide the opportunity to persistently observe the evolution of water properties, interpretation of their measurements is dependent upon knowledge of vehicle location in space and time. Without external navigational information (provided by, e.g., acoustic ranging) knowledge of vehicle location is based on vehicle hydrodynamics expressed by a flight model together with surface navigation typically provided by global positioning system (GPS) fixes. Numerous multimonth deployments of Seaglider, Spray, and Slocum vehicles (Todd et al. 2011; Cole and Rudnick 2012; Pelland et al. 2013; Timmermans and Winsor 2013) have helped improve knowledge of flight characteristics, but also highlight glider versatility and implications of piloting choices such as sampling frequency, glider vertical velocity, and glide slope angle. Sampling choices are often made to extend mission duration, but also reflect a trade-off between desired spatial and temporal resolution of ocean phenomena being observed.

Gliders nominally profile with a vertical to horizontal glide slope ratio less than unity at vertical speeds of 0.05–0.5 m s^{−1}. The effects of sampling geometry and relatively slow speed on slant profile accuracy, as compared to profiles collected by other platforms, were identified by Rudnick and Cole (2011), specifically the Doppler shifting and aliasing of higher-frequency temperature and salinity variability in measurements along depth surfaces onto large-scale structure. A shift from depth to isopycnal coordinates removes this effect and permits observations of temperature and salinity variance along isopycnals with horizontal resolution dependent on maximum profile depth and glide slope. Estimates of dynamical fields such as geostrophic velocity and relative vorticity, require consideration of horizontal density gradients along isobars and reference velocity estimates. Two examples highlight relative extremes of glider sampling strategies. Piloting choices made by Pelland et al. (2013) illustrate the approach of attempting to minimize energy use to collect a near-continuous multiyear record of the mean across-slope density and along-slope velocity structure in the northern California Current System. Seagliders repeatedly profiled to 1000 m every ~9 h along two ~200 km transects off of the Washington coast. Interpolation of the cross-slope density field and glider-inferred depth-average currents resulted in monthly snapshots of geostrophic velocities in the upper 1000 m comprising the California Current System. Aliasing of high-frequency vertical isopycnal displacements due to slow sampling speeds was reduced using Gauss–Markov interpolation along depth surfaces. Alternatively, with the approach of attempting to maximize horizontal resolution in the surface ocean, Timmermans and Winsor (2013) piloted Slocum gliders in shallow waters diving to 40 m and back every ~15 min. This sampling strategy permitted observation of both high-wavenumber along-isopycnal temperature and salinity variance, and internal wave displacements with time scales of hours. These examples illustrate different glider sampling strategies to observe large-scale density structure and finescale thermohaline variability, respectively.

More recent deployments reveal an increasingly diverse set of glider sampling strategies and objectives. Alvarez and Mourre (2012) consider optimal sampling strategies, not as a function of depth as discussed above, but instead with the aim of best resolving the horizontal structure of an evolving mesoscale field using a stationary mooring as a local reference. Repeat transects, carried out in the Mediterranean Sea, Sea of Japan, and eastern Pacific likewise demonstrate different settings in which gliders are used to estimate mixing and transport (Cotroneo et al. 2019; Wagawa et al. 2020; Jakoboski et al. 2020). These observations are all referenced to motivate consideration of the dependence of aliasing and bias on sampling pattern geometry and speed. While aliasing of high-frequency motions is unavoidable when considering horizontal gradients, piloting choices based on scientific objectives can be made to reduce it.

Deepglider, a recently developed buoyancy-driven autonomous underwater vehicle operationally similar to Seaglider but with a maximum operating depth of 6000 m, extends glider observational range considerably, prompting further consideration of the effects of sampling slowly along slanted paths in ocean. Because over 98% of the World Ocean is shallower than its maximum operating depth, Deepgliders are capable of collecting full-depth temperature and salinity slant profiles daily with vertical resolution of order meters. With the exception of wire-walking instruments deployed on moorings at multiple depths, vertical resolution of temperature and salinity measurements made by gliders is far finer than can affordably be attained by traditional moored instruments (Wunsch 1997; de La Lama et al. 2016). Despite limitations of coarse vertical resolution and full-depth deployments at relatively select locations, mooring time series of water properties have provided a description of complete eddy vertical structure limited to the gravest few vertical modes (Wunsch 1997). Vertical modes are a dynamical set of orthogonal eigenfunctions, each a time independent solution describing the vertical structures of quasigeostrophic eddies. These modes are a natural basis set often used in partitioning vertical structure into barotropic and baroclinic components in a dynamically relevant manner. Estimation of the partitioning of energy across modes, akin to a spectral analysis identifying spatial or temporal scales associated with dominant fractions of signal variance, permits consideration of energy containing scales and transfers in the vertical. The ability of Deepgliders to occupy remote regions of the ocean for durations as long as a year permits collection of more highly resolved density structure and motivates renewed consideration of complete eddy vertical structure as well as the accuracy of derived geostrophic velocity profiles based on slow sampling along slanted vertical–horizontal paths in the ocean.

High-resolution ocean numerical model simulations provide opportunity to simulate idealized glider flight across a range of speeds and glide slopes. The simulated collection of measurements permits an analysis of glider profiling efficacy analogous to that made by Rudnick and Cole (2011). We are interested in the accuracy of 1) geostrophic velocity profile estimates referenced to a glider-inferred depth-average current, and 2) isopycnal vertical displacement profile estimates, both associated with quasigeostrophic eddies. Using output from two independent ocean numerical model simulations, one in the eastern North Pacific and the other in the North Atlantic, we simulate glider sampling along zonal and meridional transects and compare glider-sampling-derived geostrophic velocity fields to instantaneous and multiday-average numerical model velocities.

This paper is organized as follows: section 2 discusses Deepglider sampling strategies employed during a 10-month Deepglider mission in the North Atlantic subtropical gyre to observe evolution of the mesoscale eddy field. Model resolution and output from the Hybrid Coordinate Ocean Model (HYCOM) and LiveOcean, a regional ocean modeling system (ROMS) model, are also detailed. Section 3 describes the simulation of glider sampling through each model, the techniques used in estimating isopycnal vertical displacement and geostrophic shear profiles, and the calculation of vertical modes. Section 4 presents results, discusses glider sampling error of velocity profiles as a function of glide slope and speed, and quantifies glider geostrophic shear estimate error along with its dependency on a ratio between eddy and high-frequency (e.g., tidal) isopycnal variability. Section 5 concludes with some advice on glider piloting strategies.

## 2. Framework

### a. Deepglider flight and observations in the North Atlantic

Deepglider is a buoyancy-driven autonomous underwater vehicle capable of profiling from the surface to 6000 m and back in approximately 1.5 days. Developed as a Seaglider variant capable of regularly profiling to the seafloor, Deepgliders employ an adapted Seaglider flight model and hydrodynamic framework originally described in Eriksen et al. (2001). Temperature and conductivity measurements are made by a SeaBird Electronics thermistor and conductivity cell aspirated by vehicle motion, processed to calculate salinity following methods described in Pelland et al. (2013). Measurements are taken at pilot specified time intervals, typically 10 s in the upper ocean increasing to 60 s deeper than 1000 m. Vertical bin averaging of measurements to a nonuniform grid, with intervals selected to obtain an approximately equal number of measurements in each depth bin, provides a nominal resolution of 1 m over the upper 150 m, 5 m between 150 and 300 m, 10 m between 300 and 1000 m, and 20 m from 1000 m to the seafloor. Using these vertical bin intervals, a typical bin-averaged temperature slant profile from the surface to 5000 m is composed of over 400 samples. Gliders also estimate a depth-average current (DAC) as the difference between GPS-tracked overground and dead-reckoned displacement. Temperature, salinity, and depth-average current accuracy are presumed 0.003°C, 0.01 psu, and 0.01 m s^{−1} (Pelland et al. 2013; Todd et al. 2011), respectively, and subsequently used in scaling simulated random instrument noise.

During descent and ascent, gliders translate horizontally with a vertical to horizontal glide slope dependent on glider pitch, buoyancy, hydrodynamic lift and drag parameters, depth, and horizontal range to a target location. Glide slope *s* is nominally 1/3, selected to roughly maximize horizontal speed for a given nominal buoyancy and vertical speed. In cross-section, successive dive–climb cycles completed along a track of constant compass heading in still water result in a sawtooth-track pattern (Fig. 1) with measurements distributed in depth and distance with time along the path. Horizontal temperature and salinity gradients along isobars or isopycnals can then be estimated from samples along the sawtooth path. The horizontal distances over which these gradients are estimated vary with depth and are a function of maximum sampling depth and glide slope. For a vertical speed of 0.075 m s^{−1} and glide slope *s* = 1/3, two full-depth slanting tracks to 5000 m, one each in descent and ascent, are completed about every 37 h and span a horizontal distance of about 30 km through the water.

Deepglider sg035, used for reference in subsequent sampling simulations, was deployed offshore Bermuda and sent to the Bermuda Atlantic Time Series station (BATS) site (31°40′N, 64°10′W) in February 2015 where it collected 350 full-depth slant profiles during a 10-month repeat survey mission (Fig. 2). Each leg of the bowtie repeat pattern was composed of two to four pairs of downward and upward slant paths, each pair labeled as a dive–climb cycle. The horizontal track of the dive–climb cycles of Fig. 1 is shown in Fig. 2, an example of a relatively straight track along a bowtie leg. The bowtie pattern was selected to permit estimation of gradients and geostrophic velocities in both the directions parallel and normal to Bermuda Rise topography. Simulations of glider sampling of ocean numerical model fields follow a similar pattern with uniformly directed transects.

### b. HYCOM

Two weeks of hourly full-depth temperature, salinity, pressure, and horizontal velocity fields were extracted along four transects, each 191 km in horizontal length, from climatologically forced and eddy-resolving HYCOM simulations near the BATS site (Fig. 3). The numerical configuration of these simulations is similar to that presented in Xu et al. (2016) and Chassignet and Xu (2017) with domain bounds near 28°S in the South Atlantic and at the Fram Strait (80°N) in the North Atlantic. In the region near BATS, the model is isopycnic with 32 vertical levels and a horizontal resolution of approximately 3.8 km (1/25°). Simulations obtained for this analysis include tides, with eight tidal constituents, and reflect on-going experiments to incorporate tides into all HYCOM simulations (details will be documented in the future).

Along the extracted slices, vertical levels are spaced at intervals of less than 10 m between the surface and 40 m, ~100 m between 50 and 1200 m, and 250 m increasing to 550 m between 1200 m and the seafloor. Model fields were interpolated to a grid with vertical spacing of 10 m from the surface to 100 m and 20 m from 100 m to the seafloor.

### c. LiveOcean ROMS

One month of hourly full-depth temperature, salinity, pressure, and horizontal velocity fields were extracted along a 200 km zonal transect from LiveOcean, a ROMS simulation in the eastern North Pacific (Fig. 4). LiveOcean (https://faculty.washington.edu/pmacc/LO/LiveOcean.html) is a regional forecast simulation with variable horizontal resolution of 0.5–1.5 km and 30 vertical levels. Ocean boundaries are forced by HYCOM ocean fields and the surface is forced by winds, solar heating, and air-sea exchange. This model includes river inflows and tides (Giddings et al. 2014). A relatively flat region of the model domain was selected with an average depth, offshore of the continental slope, of 2800 m. Along this transect, vertical levels are spaced at a variable interval of 6–175 m and were also interpolated to a grid with vertical spacing of 10 m from the surface to 100 m and 20 m from 100 m to the seafloor.

## 3. Methods

### a. Simulated glider tracks

The slanting paths of simulated steady glider flight along model zonal and meridional transects (Figs. 3, 4) were generated by selecting a starting time and location along a model transect, a fixed vertical speed, maximum dive–cycle depth, and glide slope. Sampling time intervals were set by the specified vertical speed and vertical grid to which fields in both models were interpolated. Beginning at the specified starting surface location, model fields of temperature and salinity were linearly interpolated at each vertical grid point to the time and horizontal location of the simulated glider path at that depth, determined by speed and glide slope specifications. Interpolation continued to the specified maximum depth and likewise to the sea surface during the climb portion of each cycle. The maximum depth attained on each profile is 10 m above the model bottom. For each simulated dive–climb cycle, a depth-average current estimate is calculated as the average of all cross-track model velocities within the spatial and temporal limits of the dive–climb cycle. These depth-average current estimates are subsequently used in referencing glider-estimated cross-track relative geostrophic velocity profiles. Sampling error is added to each measurement of temperature, salinity, and the depth-average current using random samples from a normal distribution with zero mean and a standard deviation equal to half of the instrument accuracy estimates described in section 2a. This procedure is repeated until the desired number of dive–climb cycles is reached or the glider reaches the end of the model domain.

Along each model transect, sets of four dive–climb cycles were simulated for starting surface positions of 10, 40, 70, and 100 km from the transect’s edge, and for starting times at the earliest time of model output and at subsequent 72-h intervals. Simulated glider profiling beginning at these starting surface positions and start times results in sampling of a variety of dynamic model features including surface intensified eddies with velocities of over 0.5 m s^{−1} and radii of over 40 km as well as the internal wave field in the absence of strong mesoscale features. Simulated glider sampling at these starting locations and start times was repeated for four different realistic glider vertical speeds, 0.06, 0.075, 0.1, and 0.2 m s^{−1} and vertical to horizontal glide slopes of *s* = 1/2 and *s* = 1/3.

Simulated glider slant profiles of temperature and salinity along the specified vertical grid are used in the calculation of neutral density *γ*^{n} (Jackett and McDougall 1997). Slant profiles of neutral density are then used to estimate isopycnal vertical displacement about time mean depths as well as horizontal density gradients along isobars.

### b. “W”/“M” sampling: Isopycnal vertical displacement and geostrophic velocity profile estimation

Accuracy of the technique employed to estimate isopycnal vertical displacement and cross-track geostrophic velocity profiles at station BATS is explored through simulation of glider sampling along transects of HYCOM and LiveOcean model output. While glider slant profiles completed along each transect leg of the BATS bowtie pattern do not strictly follow a line in plan view (Figs. 2, 5), the technique used to estimate vertical displacement and a horizontal density gradients as a function of depth and along linear paths is used in practice on sets of dive–climb cycles completed while glider heading does not deviate by more than 40° from an average heading.

In this framework, each displacement and velocity profile estimate is made using samples from two dive–climb cycles that form a pattern similar to the letter “W” or similarly from a pair of consecutive climb–dive cycles forming a pattern similar to the letter “M” (Fig. 6). Displacement profile estimates are made using four samples at each depth and geostrophic velocity profile estimates can be derived following the estimation of horizontal density gradients using thermal wind. Each horizontal gradient estimate is calculated from four samples at the same depth along respective slant paths of the pairs of dive–climb (W) or climb–dive (M) cycles. The horizontal distance over which gradients are estimated is double at the surface what it is at the bottom for W sampling and conversely for M sampling. The estimated vertical profiles of displacement and along-track horizontal density gradients are located at the center of each pattern, the junctions between the pairs of cycles used.

Because actual glider transects rarely are precisely straight in plan view, horizontal density gradients at each depth are estimated applying a least squares two-dimensional planar fit to four density measurements. The two-dimensional horizontal gradient is then projected onto the mean along-track direction (Fig. 5). For simplicity, simulated glider fight paths in LiveOcean and HYCOM are along zonal or meridional transects and errors introduced as a result of deviations in heading are not considered here.

A single simulated glider transect composed of four dive–climb cycles permits estimation of five vertical isopycnal displacement and absolute geostrophic velocity profiles located at the junctions of all dive–climb and climb–dive pairs of slant profiles. Each displacement and velocity profile is computed from two successive dive–climb (climb–dive) simulated tracks to form a vertical profile estimate centered at the midpoint of the sampling W (M) window. (Fig. 6). At grid depths *z*_{i}, average vertical displacement *ξ* in meters, across two consecutive dive–climb or climb–dive cycles, is defined as

where *γ*^{n}(*z*) is an average density profile of the four contributing slant profiles and $\gamma n\xaf\u2061(z)$ is a temporal average of all model density profiles at the midpoint location of the W or M pattern along the transect. This temporal average of model density profiles spans the entire record of model output obtained for this analysis, a 2-week period in HYCOM and one month in LiveOcean, and is meant to reflect the mean density field atop of which eddy perturbations occur.

Because simulated glider sampling is along a constant heading, one-dimensional along-track horizontal density gradients ∂*γ*^{n}/∂*x* are calculated at each depth, where *x* is along-transect distance. While these horizontal density gradients reflect the combined effect of all model dynamics, we seek to capture quasigeostrophic density gradients, presumed evolving on time scales greater than the approximate time required to complete two dive–climb or climb–dive cycles and unchanging throughout this period over a horizontal distance equal to the greatest distance between contributing measurements. The assumption that these gradients are entirely quasigeostrophic is incorrect, and aliasing of higher-frequency and wavenumber dynamics onto these estimates is subsequently explored, but it permits the estimation of cross-track geostrophic velocity shear as

where *υ*(*z*) is cross-track geostrophic velocity and is positive to the left of *x*, *g* is gravitational acceleration, *ρ*_{0} a reference density, and *f* the local Coriolis parameter. Following estimation of horizontal density gradients at each depth for each dive–climb (climb–dive) cycle pair, the vertical shear of the geostrophic velocity is calculated, integrated from the seafloor to the surface, and referenced to the glider-estimated depth-average current to obtain a full-depth profile of absolute geostrophic velocity. The location along the transect of each profile is set at the midpoint location of contributing density profiles.

In what follows, the terms glider velocity profile and glider displacement profile are taken to mean estimates of vertical profiles of these quantities estimated from pairs of consecutive simulated glider cycles (dive–climb for “W” sampling and climb–dive for “M” sampling).

### c. Vertical wavenumber spectra

Vertical structures of glider and model isopycnal vertical displacement and horizontal velocity profiles are explored through the use of normal vertical modes. Quasigeostrophic dynamics, descriptive of the mesoscale eddy field, permit a separation of variables solution for horizontal velocities *u*, *υ*(*x*, *y*, *z*, *t*) in which eddy vertical structure of slowly evolving geostrophic velocity can be described by a set of *m* normal modes *ϕ*_{m} satisfying

where *N*^{2}(*z*) is the Brunt–Väisälä frequency calculated using an average background density profile, *λ*_{m} the *m*th eigenvalue solution (s m^{−1}), and *ϕ*_{m} the *m*th eigenfunction describing horizontal velocity structure (Wunsch 1997). Numerical solutions to Eq. (3) are obtained for the first 30 modes with flat bottom and free surface boundary conditions (the first three of which are depicted in Fig. 7). These functions, the vertical modes, are normalized per unit energy such that

where *H* is the full depth of the water column. This normalization permits water-column average kinetic energy to be expressed solely as a function of mode amplitude. Vertical modes are first projected onto *i* individual cross-track glider and model horizontal velocity profiles generally expressed as

to obtain mode amplitudes *α*_{mi}. These profiles are squared and vertically averaged:

By the orthogonality of the modes, cross terms in the product of sums drop out giving

The result is water-column average kinetic energy ⟨*υ*^{2}⟩, across the sampling domain, expressed as a function of mode number and amplitude ⟨*α*_{m}⟩:

Due to the geometry of extracted model transects and the technique used to compute cross-track geostrophic velocity, this kinetic energy estimate includes only a single component of horizontal velocity, its cross-track value. Anisotropy of the velocity field across the model domain and in time is not considered.

Estimates of isopycnal vertical displacement are similarly used to explore the partition of potential energy across vertical modes. Isopycnal vertical displacements, presumed to largely reflect mesoscale eddy isopycnal perturbations about a mean background state, are defined as

where $\beta m$ mode amplitudes are calculated for *i* displacement profiles. Actual and simulated glider, as well as model, isopycnal vertical displacement profiles capture vertical motions caused by some combination of internal wave, tide, and eddy displacements. The aliasing of higher frequency and wavenumber internal tide displacements onto eddy structure will be subsequently addressed. The relationship between the vertical structure of geostrophic velocity $\varphi m\u2061(z)$ and isopycnal vertical displacement $1/N2\u2061(z)\u2061(\u2202\varphi m\u2061(z)/\u2202z)$ is derived following Wunsch and Stammer (1997). Water-column average potential energy, as a function of mode number, is then

having substituted into Eq. (3) and normalized per unit energy. Use of basis functions satisfying Eq. (3) to describe vertical structure of velocity and displacement permits an analysis of variability in the partitioning of average water-column energy across modes. Profiles and modal projections can be grouped and averaged to compare energy partitioning across modes at different locations, at different times, for different glide slopes, or for different glider vertical speeds. Individual spectra, grouped by glide slope and/or glider vertical speed, can also be averaged and compared to reveal displacement and velocity mode structure resolution associated with different glider W/M sampling patterns.

## 4. Results: Glider–model comparison and discussion

Glider displacement and velocity profile estimates are paired with numerical model displacement and velocity profiles to determine the effect of W/M sampling, glide slope, and glider vertical speed on displacement and velocity profile estimate accuracy (Fig. 7). Model displacement profiles *ξ*_{model}(*z*) are calculated using Eq. (1) with each model density profile $\gamma modeln\u2061(z)lw,tw\xaf$ [or $\gamma modeln\u2061(z)lm,tm\xaf$] defined as the horizontal and temporal average of instantaneous model density profiles $\gamma modeln\u2061(x,\u2009z,\u2009t)$ at horizontal grid points and times occurring between the first and last samples of each W (or M) pattern (Fig. 8). A typical along-transect distance *l*_{w} (or *l*_{m}) is *O*(50) km and duration *t*_{w} (or *t*_{m}) is *O*(5) days. To obtain a model density anomaly, the numerator in Eq. (1), this profile is differenced from an overall background model density profile $\gamma modeln\u2061(x,\u2009z)t\xaf$ defined as the temporal average of all instantaneous model density profiles in the time series (approximately one month in ROMS and two weeks in HYCOM) at the horizontal location *x* of the *W* or *M* profile. Model cross-track velocity profiles $\upsilon model\u2061(z)\xaf$ are defined as the horizontal and temporal average of all cross-track velocity profiles within each *l*_{w} (or *l*_{m}) and *t*_{w} (or *t*_{m}). The horizontal and temporal averaging of these model profiles permits an appropriate comparison of model output to glider-derived profiles, each calculated from samples taken across tens of kilometers and over multiple days. This averaging also serves as a low-pass filter, reducing high-wavenumber and high-frequency motions, as does the W/M method itself. Glider–model pairs of displacement and velocity profiles are differenced to obtain average and average square error as a function of depth, glider vertical speed, and glide slope.

The average difference between 364 glider sample and LiveOcean average model velocity profile pairs is less than ± 0.02 m s^{−1} with maximum absolute differences of ~0.05 m s^{−1}. These differences are not a strong function of glider vertical speed (Fig. 9). The standard deviation of this glider–model error is enhanced in the upper 250 m, with mean square error between glider and model velocities multiple times greater in the upper 250 m than between 250 and 2500 m (Fig. 10). This enhanced upper ocean error can be attributed to the difference between *t*_{w} (or *t*_{m}) and the advective time scale of near-surface velocities. Model near-surface velocities evolve on time scales much shorter than the W/M sampling period. Additionally, glider finite difference estimates of horizontal density gradients in the upper 250 m are calculated from samples either spanning a greater horizontal distance than estimates made at middepths or from samples in two sets, spaced closely together in each (for respective W and M patterns). A modest increase in the standard deviation of glider–model velocity error in the deepest 500 m of the water column is similarly the result of sampling pattern geometry with samples spaced either at greater horizontal distances apart or in two sets of closely spaced samples. Neither glider vertical speed nor glide slope has a significant effect on glider–model velocity error, but these errors are slightly reduced for slower vertical speeds and smaller horizontal glide slopes (*s* = 1/3 vs 1/2). Sampling simulations were rerun for a larger depth-average current error of 0.02 m s^{−1}, the upper bound of estimated error from Rudnick et al. (2018), with minimal effect on mean velocity profile error.

Similar patterns of average glider–model velocity profile error are observed in HYCOM simulations of glider sampling. While HYCOM average glider–model velocity profile error remains less than 0.05 m s^{−1} throughout the water column, error standard deviation in the upper 1000 m is roughly 10 times larger than values deeper than 1000 m (Fig. 11). This marked increase in error variability can be attributed to highly variable upper-ocean along-transect velocities of order 0.4 m s^{−1}, absent in LiveOcean simulations, and enhanced internal tide isopycnal vertical displacements with time scales shorter than the W/M sampling period. Additionally, maximum sampling depths in HYCOM simulations are ~60% greater than those in LiveOcean simulations. This results in an increased W/M sampling period, with samples in the upper 1000 m having a greater time interval between initial and final simulated glider samples. These two factors, highly variable and large-magnitude along-track velocities, as well as a longer W/M sampling period, result in greater glider–model velocity error variability. Both average error and error standard deviation are reduced in the upper 500 m by approximately 50% for glider vertical speeds of 0.2 m s^{−1} (Fig. 11) and overall for a glide slope *s* = 1/3 (Fig. 12).

Displacement and velocity modes were projected onto glider and model profile pairs to obtain mode amplitudes and water-column average potential and kinetic energy as a function of mode number. Energy spectra were calculated from profiles associated with each glider vertical speed and glide slope. Agreement between glider and model water-column average potential and kinetic energy spectra, partitioned by vertical mode, is reflective of the vertical spectral resolution permitted by the W/M sampling framework for various vertical speed and glide slope combinations. Across both models, all speeds, and two glide slopes, glider potential energy spectra reproduce model spectra from the first through thirtieth baroclinic mode (Figs. 13a,c, 14a,c). The local horizontal and temporal averaging of model density profiles across *l*_{w} (or *l*_{m}) and *t*_{w} (or *t*_{m}), permits a relevant comparison of glider and model displacement profiles, reducing model potential energy across all mode numbers while retaining the same spectral shape as the model average potential energy spectrum calculated from instantaneous density profiles. This is similarly the case for average and instantaneous model velocity profiles used in calculating kinetic energy spectra.

In LiveOcean simulations, glider kinetic energy spectra do not significantly vary with glider vertical speed (Figs. 13b,d). For both glide slopes, glider and model spectra are indistinguishable from the barotropic mode through the ~eighth baroclinic mode. Spectra derived from slant profiles with a glide slope of 1/2 extend this agreement through mode ~13. This suggests that the horizontal distance spanning a W/M pattern, across which contributing measurements are separated and dependent on glide slope and water depth, limits vertical resolution at modes higher than 8. This occurs despite the fact that horizontal scales (respective Rossby radii of deformation) associated with modes 2 through 8 are shorter than the distance spanned by each W/M pattern.

In HYCOM simulations, similar patterns are observed, with the deviation between model and glider kinetic energy spectra occurring at mode 11 for a glide slope *s* = 1/3 and mode 20 for a glide slope *s* = 1/2 (Figs. 14b,d). Additionally, increased glider vertical speeds result in better agreement between glider and model kinetic energy levels at low modes. Slow glider vertical speeds result in velocity profiles with speeds greater than average model speeds and a more energetic kinetic energy spectrum. With both models limited to a vertical resolution of 30 vertical levels, expectation of the highest resolvable mode structure should be limited to at most mode 30.

While the projection of dynamical vertical modes onto glider-derived displacement and velocity profiles reveals the highest-mode vertical structure that glider sampling can resolve, glider–model velocity differences can be explicitly quantified considering the aliasing of high-frequency isopycnal variability onto glider-derived horizontal density gradient estimates. Both model simulations were selected because they contain internal tide dynamics and offer temporal resolution much greater than an inertial period. These selection criteria ensure the aliasing of high-frequency motions onto glider slant profiles, each completed over a multihour period. Comparison of model energetics reveal LiveOcean to contain significantly greater high-frequency, high-wavenumber energy than HYCOM. Frequency and zonal wavenumber spectra at 1000 m reveal LiveOcean to be more energetic than HYCOM at time scales between ~2 and 24 h and across all zonal wavenumbers (Fig. 15). Despite these differences, we expect simulated glider sampling in both models, especially LiveOcean, to readily alias motions at high frequencies and across small horizontal scales.

Over the course of an approximately 30-h W/M sampling period in LiveOcean, isopycnals are vertically deflected at each horizontal grid point on hourly time scales with frequencies shorter than a 24-h period, reflective of the internal tide (Fig. 16). These deflections about a multiday-mean isopycnal vertical structure are aliased into glider horizontal density gradient estimates. The severity of this aliasing is a function of both internal tide amplitude and mean horizontal density gradient magnitude. Mean horizontal density gradients at each depth, unchanging over a multiday period, reflect mesoscale features in which geostrophic balance is dominant. Considering these gradients a signal, and higher-frequency isopycnal vertical displacements a noise, a model signal-to-noise ratio is estimated for each glider–model profile pair at each depth. At the horizontal location of each glider–model velocity profile pair, and at each depth *i*, the eddy signal *g*(*z*_{i}) is defined as

where $\gamma modeln\u2061(xj,\u2009zi)tw\xaf$ [or $\gamma modeln\u2061(xj,\u2009zi)tm\xaf$] is the time average of instantaneous model density $\gamma modeln\u2061(x,\u2009z,\u2009t)$ across the interval *t*_{w} (or *t*_{m}), $\gamma modeln\u2061(zi)lw,tw\xaf$ [or $\gamma modeln\u2061(zi)lm,tm\xaf$] is the mean model density profile across *l*_{w} (or *l*_{m}) and *t*_{w} (or *t*_{m}), and *L* is the number of horizontal grid points with indices *j* that fall within *l*_{w} (or *l*_{m}). A larger value corresponds to a greater horizontal density gradient spanning the W (or M) pattern. The noise *f*(*z*_{i}) is defined as

where *T* is the number of time points with indices *k* that fall within *t*_{w} (or *t*_{m}). A higher value corresponds to larger-amplitude isopycnal vertical displacements with time scales shorter than *t*_{w} (or *t*_{m}). The ratio *g*(*z*)/*f*(*z*) then provides an estimate of the relative strength of the mesoscale eddy gradients as compared to higher-frequency isopycnal fluctuations. This quantity can then be used, with some predictive power, to test the ability of the W/M glider velocity profile estimation tool to accurately reproduce model multiday-mean velocities.

Signal-to-noise ratio estimates derived from model density profiles are used to select glider velocity estimates that should more accurately reproduce model mesoscale velocities. This hypothesis is tested on glider and model estimates of vertical shear of geostrophic velocity. For each glider–model profile pair the difference between glider and model vertical shear, estimated using Eq. (2) from linear horizontal density gradients, is calculated as percent glider–model shear error (Fig. 17). Despite the high occurrence of individual shear errors greater than 100%, vertical integration and reference to a depth-average current limits average glider–model velocity profile error to less than ~± 0.02 m s^{−1} across a majority of the water column.

The mean of glider–model shear error estimates with a signal-to-noise ratio greater than unity shows a reduced error of ~30% or less at nearly all depths from a total average of over 100%. Reduction in average glider–model shear error for measurements associated with a signal-to-noise ratio greater than unity confirms that glider velocity profile accuracy is improved when sampling mesoscale features with strong horizontal density gradients and relatively weak internal tides. A priori knowledge of internal tide amplitudes and the strength of mesoscale gradients can thus aid pilot decision making in seeking to minimize the aliasing of higher-frequency isopycnal fluctuations onto mesoscale density structure. While internal tide and other high-frequency aliasing cannot be avoided, geographic surveys identifying regions of increased high-frequency energy, like Zhao (2016), can inform to some extent on expected errors. Agreement between glider and model vertical energy spectra is thus improved considering the relative dominance of internal tides or mesoscale features. This is apparent considering simulated sampling of an mesoscale eddy in LiveOcean. Glider–model kinetic energy spectra agree through higher mode number when considering velocity profiles sampling only a mesoscale eddy (Fig. 18). For a glide slope *s* = 1/2, LiveOcean and glider kinetic energy spectra agree through mode 30.

## 5. Conclusions

Results of glider sampling simulations across two models, four glider vertical speeds, and two glide slopes suggest that the W/M framework can accurately resolve quasigeostrophic eddy isopycnal vertical displacements, horizontal density gradients, and geostrophic velocities. In LiveOcean, this accuracy is primarily a function of mesoscale eddy strength and glide slope, while in HYCOM, faster glider vertical speeds decrease glider–model velocity error in the upper 500 m and correspond to better agreement between glider and model kinetic energy spectra. These results suggest that as maximum glider sampling depth increases, near-surface velocity error increases, but in proportion to the magnitude of near surface advective velocities. As Deepglider missions increase in number, consideration of these factors can help in achieving mission goals. Glider sampling strategies always reflect a trade-off between energy usage and desired spatiotemporal resolution. These simulations quantify the effect of some piloting choices and show that increased glider vertical speeds are not always necessary to resolve quasigeostrophic features that evolve on multiday or longer time scales.

Across all simulations, vertical structure of eddy density anomalies and geostrophic velocities, from the barotropic through at minimum the eighth baroclinic mode, is accurately reproduced by glider slant profiling and W/M velocity profile estimation. Furthermore, the partitioning of energy across these vertical modes, containing over 95% of observed eddy mechanical energy, reveals a common spectral pattern in which energy decreases with increasing mode number *m*, or scaled vertical wavenumber, following a log–log linear slope proportional to m^{−3}. This slope agrees with geostrophic turbulence predictions on the partition of energy across modes within the enstrophy inertial range and suggests that model fields reflect quasigeostrophic dynamics.

These simulations of glider-slant profiling and geostrophic velocity profile error analyses lend confidence to the application of this framework to full-depth Deepglider observations, detailing the evolution of quasigeostrophic motions with greatly increased vertical and temporal resolution. Analyses of these new observations using this framework thus provide new opportunity to test theoretical predictions of vertical energy partitioning and transfer, as well as identify relationships to the distribution of mesoscale eddy kinetic energy across horizontal scales, of which observations are relatively more numerous.

## Acknowledgments

This work was funded by the National Science Foundation (Award 1736217). We thank Xiaobiao Xu and Eric Chassignet at Florida State University for providing accesses to HYCOM output from ongoing experiments that include tidal forcings as well as for saving four-dimensional output at hourly time steps. Their advice and assistance helped ensure proper use of the relevant model fields. We would also like to thank Parker MacCready at the University of Washington for providing access to LiveOcean ROMS output and assistance in interpretation of coastal and model boundary dynamics.

## REFERENCES

*J. Atmos. Oceanic Technol.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*Earth Syst. Sci. Data*

*Dyn. Stat. Climate Syst.*

*IEEE J. Oceanic Eng.*

*J. Geophys. Res. Oceans*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Geophys. Res. Oceans*

*J. Atmos. Oceanic Technol.*

*J. Atmos. Oceanic Technol.*

*Cont. Shelf Res.*

*J. Geophys. Res.*

*J. Mar. Syst.*

*J. Phys. Oceanogr.*

*Rev. Geophys.*

*J. Climate*

_{2}internal tides

*J. Phys. Oceanogr.*