Volumetric measurements of wake impulse and kinetic energy for evaluating swimming performance

Volumetric flow measurements are a valuable tool for studies of aquatic locomotion. In addition to visualizing complex propulsive behaviors (e.g., highly three-dimensional kinematics or multi-propulsor interactions), volumetric wake measurements can enable direct calculation of metrics for locomotive performance including the hydrodynamic impulse and wake kinetic energy. These metrics are commonly used in PIV and PTV studies of swimming organisms, but derivations from planar data often rely on simplifying assumptions about the wake (e.g., geometry, orientation, or interactions). This study characterizes errors in deriving wake impulse and kinetic energy directly from volumetric data in relation to experimental parameters including the level of noise, the flow feature resolution, processing parameters, and the calculation domain. We consider three vortex ring-like test cases: a synthetic spherical vortex with exact solutions for its impulse and energy, volumetric PIV measurements of a turbulent vortex ring, and volumetric PIV measurements of a turning fish. We find that direct calculations of hydrodynamic impulse are robust when derived from a volumetric experiment. We also show that kinetic energy estimates are feasible at experiment resolutions, but are more sensitive to experiment design and processing parameters, which may limit efficiency estimates or comparisons between studies or organisms.


Introduction
Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) are frequently used experimental techniques for studying aquatic locomotion by directly measuring velocity fields. In addition to flow field visualization, PIV and PTV measurements can be used to derive metrics for momentum and energy transfer from an organism to the surrounding fluid (e.g., Müller et al. 1997;Drucker and Lauder 1999;Stamhuis et al. 2002;Epps and Techet 2007). Metrics including the hydrodynamic impulse and wake kinetic energy can be used to make comparisons between forward swimming and turning (Wolfgang et al. 1999), at varying surrounding flow speeds (Drucker and Lauder 1999), or between organisms (e.g., Drucker and Lauder 2000;Tytell 2007;Dabiri et al. 2010). Wake metrics can also be used to evaluate how kinematic variations affect swimming performance (Akanyeti et al. 2017) and the contributions of multiple propulsors (Tytell 2006).
Volumetric flow measurement techniques such as defocusing digital PIV/PTV (DDPTV, Pereira et al. 2000), tomographic PIV (TomoPIV, Elsinga et al. 2006), synthetic aperture PIV (SAPIV, Belden et al. 2010), and Lagrangian particle tracking (Shake-the-Box or LPT, Schanz et al. 2016), are a transformative tool for studies of biological and bioinspired locomotion through fluids. For instance, in studies of fish swimming, volumetric measurements can reduce constraints on animal behavior during experiments (Adhikari and Longmire 2012a) or enable simultaneous imaging of multiple fins (Mendelson and Techet 2020).
Volumetric measurements can also relax assumptions about the geometry and orientation of wake features. When live animals are studied using planar velocimetry, it is unlikely that hydrodynamic features will align perfectly with the measurement plane, but alignment and symmetry are commonly assumed when deriving quantities such as the hydrodynamic impulse (e.g., Epps and Techet 2007). For a volumetric PIV experiment, instead of assuming alignment, Mendelson and Techet (2015b) highlight multiple reference frames for wake impulse analysis that can be selected based on organism kinematics or measured vorticity vectors. Additionally, with planar data, quantities such as wake kinetic energy may be determined per depth (e.g., Dabiri et al. 2010), assuming a consistent measurement location relative to the organism, instead of for the complete flow field. As an alternative to planar derivations, quantities such as the impulse and kinetic energy can be calculated directly from volumetric velocity and vorticity data.
Deriving wake metrics directly from volumetric velocity or vorticity data has several potential error sources and experiment design trade-offs. As with all quantities derived from experimental velocity measurements, noise and error propagation are a significant consideration. Resolution may also have significant influence, especially in cases where spatial resolution must be balanced with achieving a sufficiently wide field of view to image a freely swimming organism or multiple propulsors. For a 3D particle reconstruction, the uncertainty in the depth dimension (hereafter referred to as the Z-direction) can also be higher than the others due to particle elongation (Discetti et al. 2013), suggesting that flow feature alignment relative to the measurement volume may have an effect on quantitative analysis. For PIV techniques, volumetric measurements may require larger interrogation windows which smooth the velocity field, though window overlap is also a factor in the final resolution (Tokgoz et al. 2012). Dense 3D particle tracking techniques that consider individual particle trajectories (Schanz et al. 2016) offer an alternative to cross-correlation and have also been used to generate wake measurements for evaluating propulsive behaviors (Rival and van Oudheusden 2017). Trackingbased techniques can improve spatial gradient resolution and aid in temporal differentiation. These measurements can be interpolated onto a Cartesian grid for deriving further quantities (e.g., Schanz et al. 2018), resulting in similar outputs between cross-correlation and tracking-based velocimetry techniques.
This study evaluates two common frameworks for quantifying aquatic locomotive performance in the context of volumetric measurements, including the effects of resolution, alignment, noise, and field of view. To ensure that the results are generalizable to any volumetric velocimetry technique, we focus on uniform grid velocity fields obtained by either PIV or interpolation of PTV results.

Quantities of interest
In this study, we focus on two common metrics for quantifying the performance of a swimming organism from measurements of the vortical wake: 1. Hydrodynamic impulse, and related quantities including vortex circulation and impulse time derivatives, which are often used in estimates of forces. 2. Wake kinetic energy, and related quantities including wake power and dissipation.
Both quantities can be directly calculated from volumetric velocity and vorticity fields in a control volume.

Hydrodynamic impulse
Hydrodynamic impulse ( ⃗ I ) is defined for a vortical flow as: where is the fluid density, ⃗ x is a position coordinate, ⃗ = ∇ × ⃗ u is a vorticity field derived from a velocity field ( ⃗ u ), and V is the fluid domain. ⃗ I represents the instantaneous impulse that would need to act over a finite domain to generate the vorticity field from rest. Impulse calculations can be used to evaluate the fit between a measured wake and theoretical models (Murphy et al. 2012), assess the relative contributions of multiple propulsive modes (Bartol et al. 2016), or compare wakes and organism trajectories (Mendelson and Techet 2020).
Uncertainty in impulse calculation using Eq. (1) can be caused by noise, incomplete measurement of the vorticity field, or losses as flow moves outside the measurement volume.
The total momentum ⃗ P in a 3D incompressible flow field is related to the hydrodynamic impulse as: where n is the outward-facing normal of the surfaces (S) both surrounding the volume (V) and of any bodies (Kang et al. 2017). In addition to the use of impulse itself as a metric for swimming performance, forces can be inferred from the change in vortex wake impulse and additional terms describing inviscid pressure and body force contributions (Dabiri 2005;Wu et al. 2007;Rival and van Oudheusden 2017;Kang et al. 2017;Limacher et al. 2020). Kang et al. (2017) highlights how impulse can be leveraged for force predictions from finite domain flow fields with discrete wake vortices. Limacher et al. (2020) provides a comparison of momentum-based and impulse-based methods for force estimation.
Simplified models of vortex strength and geometry are often used to determine the impulse, especially with 2D measurements. For instance, the impulse for a circular thin-cored vortex ring, a common fish wake model, can be simplified to: where R is the ring radius and is the vortex circulation. The impulse vector acts along the ring's axis, which is assumed to be in-plane for a 2D PIV measurement. Equation 3 can also be modified to account for a finite vortex core thickness (Epps and Techet 2007). (1) The circulation is defined from either velocity ( ⃗ u ) or vorticity ( ⃗ ) measurements as: where d ⃗ l is a closed contour surrounding the vortex core and dA is the contour's enclosed area. Mendelson and Techet (2015b) estimated an uncertainly of 13% for circulation and diameter-based estimates of impulse for a vortex ring measured using SAPIV. In addition to impulse calculations, the circulation itself is also often correlated with force production to assess performance (e.g., Bomphrey 2012; Imamura and Matsuuchi 2013;Langley et al. 2014). Gutierrez et al. (2016) provides an extensive discussion of the sensitivity of circulation measurements to wake topology and choice of integration domain. In addition to these considerations, the vortex ring model for impulse requires geometric simplifications that may be unnecessary with full volumetric data. As a result, we focus on direct calculation of impulse (Eq. 1) in this study.

Wake kinetic energy
Locomotive performance can also be assessed by quantifying energy lost to wake formation. The kinetic energy in a control volume can be calculated from the velocity field ( ⃗ u , with any freestream velocity removed) as: Variations to Eq. (5) can account for flows with induced velocities extending beyond the measurement volume; these formulations are often used in studies of animal flight (Hedh et al. 2020;Johansson and Henningsson 2021). Temporal derivatives of energy can be used to evaluate power. The average wake power, commonly required for efficiency metrics, can be estimated as the change in wake kinetic energy divided by stroke time (Bartol et al. 2016;Hedh et al. 2020). When propulsion has ceased or the organism leaves the measurement volume, temporal decays in kinetic energy correspond to dissipation and fluxes of kinetic energy out of the measurement volume. The viscous dissipation rate ( ) can be calculated from velocity spatial gradients: where u, v, and w are the measured velocity components and is the fluid's dynamic viscosity. The accuracy of a dissipation calculation depends on both processing parameters and the flow's physical attributes. Tokgoz et al. (2012) evaluates 3D dissipation measurements for laminar and turbulent Taylor-Couette flows, noting that increasing processing window overlap improves dissipation calculations by providing better estimates of local spatial gradients and that dissipation measurements are less accurate for more turbulent cases. Table 1 summarizes studies using volumetric velocimetry methods to evaluate aquatic locomotive performance and the quantities derived in each work. Studies that use volumetric techniques for visualization or velocity and vorticity measurements only are excluded. Organism sizes range from plankton to humans. Spatial resolution is normalized by the minimum number of vectors per organism length.

Applications
Based on their frequency of use, we evaluate deriving impulse, kinetic energy, and related quantities from a volumetric vortex-dominated flow field and characterize the effects of resolution, data processing and smoothing, field of view, and flow feature alignment with the measurement volume. We utilize three test cases, all vortex ringlike structures similar to the wakes of aquatic organisms: a synthetic vortex with exact expressions for impulse and kinetic energy, time-resolved SAPIV data of a turbulent vortex ring, and time-resolved SAPIV data of a turning giant danio fish.

Synthetic flow field
To directly assess impulse and kinetic energy errors, we employ a synthetic field: Hill's vortex (Fig. 1a), which consists of axisymmetric vortical flow within a sphere and irrotational exterior flow (Hill 1894;Akhmetov 2009). We construct the velocity field in the absolute frame, where the uniform external freestream velocity ⃗ u 0 = −u 0ẑ is removed (Fig. 1b), with u 0 ≥ 0 , so that u → 0 as r → ∞ . In this frame, ⃗ u 0 refers to the uniform translation of the vortex along the positive Z-axis. Exact expressions for the impulse and kinetic energy of the vortical region are: where R is the radius of the spherical vortex (Hill 1894). We consider the kinetic energy of just the internal vortical region (Eq. 8), since it better resembles vortex ring wakes. The external flow also extends to infinity, so the true kinetic energy for the external region would always exceed the kinetic energy calculated from a finite velocity field. This domain choice has no influence on the impulse calculation, as the external region is irrotational and does not contribute to vortical impulse. We compare the computed and (7) ⃗ I = 2 R 3 u 0ẑ , theoretical values for vortical impulse and kinetic energy for the following sources of error. Noise and post-processing. The effect of random noise in the velocity field on measurements of impulse and kinetic energy is considered by adding a corresponding field of random velocity noise, i.e. ⃗ u ← ⃗ u + ⃗ u (Fig. 1c). The noise field ⃗ u has independently generated random vectors at each grid position and in each component. Specifying an upper bound for velocity noise magnitude u > 0 , the velocity noise is uniformly random in each dimension, with range . We specify u as a proportion of u 0 . For all cases including noise, 20 random iterations are performed at each noise level.
Noise models with independently generated vectors over the grid, with uniform or Gaussian distributions, have been used in studies of lift and drag calculations, including differential and integral formulations (David et al. 2009), of pressure calculation (Charonko et al. 2010;Nie et al. 2022), and in earlier studies of vorticity computation (Foucaut and Stanislas 2002). Models with local noise scaling and/or local spatial correlation have also been commonly employed (Azijli and Dwight 2015; McClure and Yarusevych 2017), including for studies with path integrals through the measurement volume, which impulse and kinetic energy do not require. For (Eqs. 1 and 5), using a two-parameter noise model did not significantly alter the main error trends observed using the uniform, single-parameter model (additional information in Supplementary Material).
To simulate aspects of PIV post-processing, we apply smoothing to the velocity field with added noise prior to calculations of impulse and kinetic energy. Two simple smoothing filters are used: box (averaging) or Gaussian, with a kernel size of 3 × 3 × 3 . This model for post-processing assumes that large global velocity outliers would be detected and replaced prior to smoothing or deriving any quantities.
Spatial resolution. We generate Hill's vortex at various feature resolutions = R∕s , where R is the radius of the vortex and s is the spacing between velocity vectors. The feature resolution describes the number of velocity vectors along a dimension per characteristic radius. All velocity fields have uniform dimensions ( n × n × n ) and uniform spacing s in all three dimensions. The vortex is centered at the origin, and spacings were chosen such that R = ms for some positive integer m, giving a field dimension of n = 2 + 1.
Windowing. We simulate PIV windowing by first creating a more resolved original field, performing sliding averaging within a uniform window of size w in each dimension, and downsampling the field by selecting the central velocity vector in each window, moving the window according to a specified overlap ratio o.
Origin and domain selection. Impulse calculation (Eq. 1) requires a choice of origin. DeVoria et al. (2014), in studying force calculations, proposed an objective choice of origin to mitigate error, defined by residual minimization for two integral identities, including Eq. 2. Similarly, we use the magnitude of the residual of Eq. 2 to characterize errors due to origin and calculation domain selection.

Experimental data
We also compare impulse and kinetic energy analysis using two sets of time-resolved synthetic aperture PIV data: a turbulent vortex ring and a turning giant danio fish. Both experiments used a camera array of nine Vision Research Miro 310 high-speed cameras and Nikon Nikkor 50 mm lenses set to f # = 8 . Fast 3D particle volume reconstruction through refractive media was performed using the homography-fit (HF) method (Bajpayee and Techet 2017). PIV interrogation was performed using multipass cross-correlation implemented in a custom 3D version of MATPIV (Sveen 2004). Table 2 summarizes experiment parameters for each case.

Turbulent vortex ring
The vortex ring was generated by a mechanical piston moving upstream of a 12.7-mm-diameter nozzle to displace a slug of fluid. The ratio of piston stroke length to orifice diameter ( L∕D o ) was 3.7, within the range of formation numbers that produces an isolated vortex ring instead of a multi-vortex jet (Gharib et al. 1998). The Reynolds number U p L∕2 = 37,000 , where U p is the piston velocity. While this Reynolds number is much higher than most aquatic wakes, the inherent velocity fluctuations are valuable for assessing post-processing and domain selection options. Synthetic aperture PIV was performed 6D o downstream of the nozzle. The resolution of the vortex ring = R∕s ≈ 5 . The propagation velocity of the ring in the downstream (Y) direction, measured from the centroid position, was constant ( ≈ 0.75 m/s, R 2 > 0.99 ) throughout the range of times considered. Thirty-six timesteps of data were analyzed. Post-processing. Multiple smoothing operations were applied to velocity field data following outlier detection (using signal-to-noise ratio and median filters) and interpolation of missing vectors: spatial smoothing using a 3 × 3 × 3 Gaussian kernel, spatial smooth using a 3 × 3 × 3 box filter, and spatiotemporal smoothing by fitting local second-order polynomials in space and time over a region of 5 × 5 × 5 vectors in space and five timesteps (Scarano and Poelma 2009). Figure 2 shows one timestep of the results from this experiment with each smoothing option applied. The degree of smoothing and noise suppression are observably higher for the box and spatiotemporal filters.
Domain selection. Six different integration domains were considered, beginning with a 30 × 45 × 37 mm bounding box surrounding the vortex ring region ( | ⃗ | > 10% of the time-averaged maximum vorticity magnitude) and expanding outward by 1-5 vectors on each surface. The origin for all impulse calculations was held fixed at the centroid of the integration domain.
Feature alignment. The turbulent vortex ring case was used to assess flow feature alignment by isolating and comparing quantities calculated from the X-Y and Z-Y planes of the flow field, which should be identical for an axisymmetric ring. The spacing of cameras used in this experiment resulted in coarser voxel spacing in the Z-direction, with z x ≈ 2.4 (Belden et al. 2010), allowing the effects of this resolution difference on derived quantities to be quantified.

Turning giant danio
The second experimental test case is data from a SAPIV experiment on a turning giant danio. The fish had body length BL = 7.5 cm and a mass of 8 g. This test case enables evaluation of impulse and kinetic energy calculation for an actual swimming wake. Figure 3 shows the turn kinematics and wake for three timesteps and two window overlaps ( o = 0.5 and o = 0.75 ). Supplementary Figure S8 shows planar slices through the wake taken through the caudal peduncle. The wake consists of smaller preparatory and larger propulsive vortices, similar to previous observations of turning and maneuvering giant danio (Wolfgang et al. 1999;Epps and Techet 2007). The flow structures in this case were less axisymmetric than the previous two test cases and in closer proximity to other wake structures and the measurement volume edges. Additionally, this case includes a moving body in the flow, which was reconstructed using the visual hull method (Adhikari and Longmire 2012b). Additional details on the processing parameters for this case are provided in Mendelson and Techet (2015a), though new post-processing was applied for the current study. Eighty-nine timesteps of data were analyzed and t = 0 s corresponds to the first frame of PIV processing.
Edge effects. Impulse and kinetic energy were calculated for a 70 × 50 × 55 mm volume behind the fish. The wake of the turn linked with previous wake structures created as the fish swam into the turn. As a result, there was vorticity on the boundaries of the analysis volume ( Z = −5 mm region of Fig. 3). We calculate fluxes to determine any leakage of Impulse fluxes have previously been used to estimate forces from wake measurements (Limacher et al. 2020). We performed forward integration in time to determine cumulative impulse or energy gains or losses: where t is the current time in the experiment and t ′ is the variable of integration. Window overlap. The turning giant danio data are also used to further evaluate processing at multiple window overlaps ( 50% and 75% ), especially when a body is present in the flow. Data post-processing for both cases consists of a signal-to-noise ratio filter on correlation peak heights, a local median filter, and spatial smoothing with a Gaussian filter. Kernel sizes for filters were 3 × 3 × 3 vectors for the 50% overlap case ( = 0.65 for Gaussian filters) and 5 × 5 × 5 vectors for the 75% overlap case ( = 1.3 for Gaussian filters), ensuring filter kernels had the same physical dimensions and comparable weights.

Spherical vortex
The synthetic Hill's vortex was used to quantify the effects of noise, resolution, windowing, and origin selection on measurements of hydrodynamic impulse and kinetic energy.

Impulse
Noise and post-processing. As impulse is a linear transformation of the velocity field (Eq. 1), by introducing noise ⃗ u , the resultant error in impulse is: The impulse error due to noise propagation is independent of the flow field and is a linear transformation of the noise. Given that ⃗ u is randomly distributed, we expect significant cancellation within this summation. A noise model that assigns independent random noise vectors to all grid positions results in the sum of a large number of independent random variables in Eq. (13), which approximates, by the central limit theorem, a normal distribution with a confined spread (Supplementary Figure S1). For a noise model simulating the local correlation of PIV measurement errors (Azijli and Dwight 2015), which effectively reduces the number of independent random variables, a normal profile is still obtained (Supplementary Figure S2). Figure 4 shows results for impulse noise propagation at fixed feature resolutions = 10 and = 20 and origin ⃗ x 0 = (0, 0, 0) . Error is normalized by the magnitude of the theoretical impulse (Eq. 7). Impulse error due to noise is equally likely to be an over-or an underestimate, so we report its magnitude. The error profiles of all three dimensions are similar, as the propagated error is independent of the flow field (Eq. 13). The nonzero error when no noise is introduced is a baseline resolution error.
The error magnitude and its variation grow with increasing noise, but even without any filtering prior to vorticity calculation, impulse calculations at = 20 (Fig. 4b) incur errors of approximately 10% or less for noise levels as high as 300% of u 0 . This rather dramatic suppression of error from noise is due to the summation of a large number of independent random variables, which results in a normal error profile with a confined spread, according to the central limit theorem. Applying a Gaussian filter, the mean error can be reduced to under approximately 5% . A Gaussian filter can similarly reduce the mean error below approximately 10% for a resolution of = 10 . Using a box filter yields a comparable trend with even lower error. Thus, due to the linearity of impulse with respect to velocity, impulse calculations are robust to noise.
Spatial resolution. Error due to the finite resolution of the measured velocity field is inevitable. We construct Hill's vortices of various resolutions, = 4 − 26 and compute the error due to resolution alone (Fig. 5a) and with the addition of noise (Fig. 5b).
The baseline resolution error was beneath 5% , which means that for Hill's vortex, a handful of uniformly sampled vectors are sufficient to estimate the vortical structure for computing impulse. Without applying a filter, the error generally decreases with higher resolution. At lower resolutions, smoothing results in an overestimate rather than an underestimate of impulse. At relatively high resolutions ( ≈ 15 in Fig. 5a), error decreases monotonically, though the improvements in both the filtered and unfiltered cases diminish with increasing resolution. The initial dip and subsequent resurgence observed at low resolutions are likely the result of imperfect sampling and unlikely to generalize to all flows.
When noise is considered by adding u = 1.5u 0 , a more monotonic trend of decreasing error magnitude with increasing resolution is observed (Fig. 5b). The mean error as well as the variation diminish as the flow features become more finely resolved and random noise more effectively suppressed by filters without neutralizing local features. For our synthetic test case, using a Gaussian filter, a moderate resolution of = 6 is sufficient to reduce the mean error beneath 10%.
Windowing effects. We compare results from different window sizes and overlap ratios using a fixed initial resolution = 160 . To avoid truncation effects, we constructed the initial field with a margin of the maximum window size (64 voxels) in the irrotational region on each edge. The underestimation of impulse due to more extensive smoothing in larger windows is shown in Fig. 6a. The downsampled fields have resolutions = 20, 10, 6.67, 5 for window sizes of 16, 32, 48, and 64 voxels, respectively, and 50% window overlap.
Overall, the error is below 10% for all window sizes considered at 50% overlap. For the 48 and 64 voxel cases, which correspond to = 6.67 and = 5 , windowing results in a greater underestimate of impulse than spatial resolution alone (Fig. 5a). When smoothing is applied to the downsampled data, it imparts some vorticity to the external irrotational region, resulting in higher impulse values and, for the 16 and 32 voxel window sizes, an overestimate instead of an underestimate.   Figure 6b shows that increasing window overlap to o = 0.75 yields the most accurate estimates of impulse across all window sizes. The resolutions of the windowed fields are = 13.3, 6.67, 4.44, 3.33 at o = 0.25 and = 40, 20, 13.33, 10 at o = 0.75 . The 25% overlap, 32 voxel window case and the 50% overlap, 48 voxel window case have the same resolution ( = 6.67 ), but the 32 voxel case has less error. Reduced window overlaps are most detrimental with larger windows, but with some variations in the amount of underestimation (e.g., the lower error for w = 64 voxels and 25% overlap than the w = 48 voxels case at the same overlap). This variation may be due to the exact placement of velocity vectors relative to the vorticity distribution, which is not something that can be easily designed for in an experiment where the wake structure is not known beforehand. These results demonstrate that increasing window overlap when larger interrogation windows are used improves estimates of impulse.
Origin selection. For a noisy measured velocity field, the choice of origin ( ⃗ x 0 ) in calculating impulse can be influential. The error contribution from each fluid element is scaled by its relative distance to the origin (Eq. 13). Assuming comparable noise across the region in which impulse is calculated, a central location which minimizes the overall distance to all fluid elements within the vortex is desirable as the origin (Supplementary Figure S3 illustrates an example). Alternately, for a specific noise profile, an error-reducing origin can be determined by minimizing the residual ( |⃗| ) of Eq. (2) We compare centroid and objective origins for Hill's vortex at = 20 with noise levels from 0% to 300% . The mean errors with objective origins, which are fitted to the specific noise profiles of each trial, are lower than the centroid origin errors, especially at significant levels of noise (Fig. 7a). However, improvements in error using an objective origin are at maximum 3 − 4%.
The momentum equation residual ( |⃗| ) and the impulse error were correlated. Best-fit lines for the objective origin and centroid data across all noise levels are shown in Fig. 7b. The correlation between the residual and the impulse error was stronger for the centroid origin cases. While impulse errors were reduced using the objective origin, there was more variation in the correspondence between the residual and the impulse error.
While the objective origin can lower error, improvements are small and comparable to other strategies for working with noisy data. For instance, if we were to perform a smoothing operation prior to calculating impulse with the natural origin, we obtain levels of error slightly lower than those from using objective origins without previous filtering of the velocity field. Optimizing for the objective origin after smoothing the velocity field still yields a marginal reduction, but on the order of less than a percent and only at high levels of noise.

Kinetic energy
Noise and post-processing. Unlike impulse, kinetic energy is not a linear function of velocity. Propagating noise in equation 5 yields The first term depends on the noise alignment with the actual velocity and is subject to cancellation within the integral, while the quadratic term produces an overestimate. At levels of noise u ≳ u , we assume that the latter term will dominate. Due to this quadratic growth, at high noise levels, postprocessing is necessary prior to kinetic energy calculations. For example, using a Hill's vortex with = 20 , the unfiltered error reaches more than 20% for a noise level of 1.0u 0 and more than 100% error for a noise level of 2.5u 0 (Fig. 8).
As with impulse, we construct Hill's vortices at resolutions = 10 and = 20 and subject them to noise levels up to 3.0u 0 . The kinetic energy of the vortical region (within the sphere) with the noisy velocity is compared to the exact solution (Eq. 8).
(14) The extent of error after filtering is very contained (Fig. 8b-c), with errors below 10% until noise exceeds 2.0u 0 . At the higher resolution, using a box filter, which more forcefully removes random noise than the Gaussian, a global bound of 6% mean error is observed for up to 3.0u 0 noise. The initial decay and subsequent rise of the mean error magnitude (Fig. 8b-c) is due to the initial underestimation of kinetic energy due to resolution and smoothing, which is then numerically compensated with the energy increase by introducing noise (the | ⃗ u| 2 term in Eq. 14). The range of errors observed across multiple trials is also smaller than impulse. Unlike impulse, which is a vector that depends on the velocity direction and the position of the fluid elements, kinetic energy is a scalar quantity which depends only on the velocity magnitude. The lack of direction and position dependence may explain the smaller variation in error as compared to impulse. Spatial resolution. We characterize resolution effects on kinetic energy from = 4 − 26 , including baseline resolution error (Fig. 9a) and the error when exposed to noise with u = 1.5 u 0 (Fig. 9b). As with impulse, a handful of velocity vectors uniformly sampled are sufficient for an accurate estimate of the kinetic energy without noise (Fig. 9a). Smoothing actions contribute to significantly higher errors at low resolutions and converge to the original error as resolution increases (Fig. 9a). For = 12 , resolution errors are confined within 5% regardless of filters.
When noise is added, the kinetic energy error exhibits an initial monotonic climb to near 0 and subsequently increases, approaching a plateau as resolution increases. At a given level of noise, there is a resolution at which the underestimation from the discretized field and smoothing is balanced by the overestimation from the noise. In an experiment, these effects are likely difficult to decouple. Windowing effects. Windowing effects in kinetic energy calculations are characterized for the same range of resolutions, window sizes, and overlaps as impulse. As a smoothing operation, windowing produces underestimates of kinetic energy (Fig. 10a). The levels of error are contained, but slightly higher than for impulse (Fig. 6). As with impulse, the smoothing effects of large window sizes cause greater errors than due to spatial resolution alone (Fig. 9). Overlap (Fig. 10b) has less effect on kinetic energy compared to impulse, likely since no gradients are taken and the window size dominates smoothing effects.

Impulse
The turbulent vortex ring case allowed us to assess impulse and kinetic energy calculations for an axisymmetric wake with multiple experimental limitations (i.e., inherent unsteadiness, limited depth and spatial resolution). Figure 11 shows each component of the impulse over time for the smallest ( 30 × 45 × 37 mm) integration volume. All impulse results are normalized by the slug impulse of the piston o LU p and is determined based on the nozzle diameter D o , stroke length L, and piston velocity U p (Gharib et al. 1998). The normalized axial impulse I * y was initially around 0.6, lower than reported by Gan and Nickels (2010) for a similar Re case measured using 2D PIV with ≈ 14 . This disagreement is likely due to the lower resolution and dynamic range of the SAPIV measurements. For all smoothing options, the volumetric impulse calculation was more consistent over time than the impulse calculated using planar data extracted from the volumetric velocity field and an axisymmetric vortex ring model (Supplementary Figure S5 and calculation procedure in Supplementary Information). The axial impulse was also observed to decay slightly over the timescale considered, which has been previously reported for turbulent vortex rings by Gan and Nickels (2010). The impulse in the X-and Z-directions was approximately zero.
Spatial smoothing operations did not have a noticeable effect on turbulent vortex ring impulse (Fig. 11). As a result, we focus on the Gaussian-smoothed and spatiotemporally smoothed data for further comparison. Figure 12a-b shows the normalized axial impulses I * y for three of the integration volumes considered using two smoothing options. Supplementary figure S6 shows the other three integration volumes, which were generally consistent with the initial (smallest) volume. With Gaussian smoothing, the larger volumes exhibited larger impulse fluctuations, likely due to additional vorticity within the domain amplified by farther proximity from the origin at the centroid. Expanding the volume did not consistently add additional impulse beyond the initial bounding box, illustrating that impulse calculation is spatially compact and can be performed where flow structures may be in close proximity to each other. Spatiotemporal smoothing Markers denote timesteps, solid lines denote Gaussian smoothing, and dotted lines denote spatiotemporal smoothing. c Kinetic energy of the X-and Z-velocity components over time for the Gaussian and spatiotemporal smoothing cases and the 30 × 45 × 37 mm volume ( Fig. 12b) reduced the range of fluctuations seen over time.
With spatiotemporal smoothing, expanding the integration volume also resulted in generally reduced axial impulse, most noticeably at t < 0.005 s. Contributions to the axial impulse from the Z-Y plane (X-vorticity components) were lower magnitude than the X-Y plane (Z-vorticity components) (Fig. 12c), possibly due to higher smoothing in the Z-direction resulting in an underestimate of gradients. For the initial integration volume (tightly surrounding the vortex), the Z-Y plane contributions were approximately 80% of the X-Y plane contributions. If we assume that the Z-Y and X-Y plane contributions to the total axial impulse should be equal, this is equivalent to a 10% ± 1% possible underestimate of total impulse due to 2.4× reduced Z-direction particle resolution and 1.2× larger interrogation windows in the Z-direction.
The extent of the discrepancy between X-Y and Z-Y planes varied with integration volume but did not vary substantially with smoothing operation. Expanding the measurement volume lowered the Z-Y plane contribution and increased the X-Y contribution, further suggesting that the integration volume should be as close as possible to the vortex boundaries to mitigate the effects of additional vorticity elsewhere in the measurement volume.
Errors due to origin selection were characterized for each integration volume using the magnitude of the residual of Eq. 2 ( |⃗| ), also normalized by the slug impulse (Fig. 12d). The mean residual across all timesteps for the four smallest integration volumes was approximately 1% of the slug impulse (for comparison, impulses measured in the experiment were 50 − 60% of the slug impulse). For smaller integration volumes, spatiotemporal smoothing reduced the range of residuals. For wider integration volumes, spatiotemporal smoothing resulted in a wider range of residuals and a higher mean residual (up to approximately 10% of the total vortex ring impulse). The lower residual for tighter integration volumes again highlights the relative spatial compactness of impulse calculations, which makes them well-suited for finite domain measurements.

Kinetic energy
Unlike the impulse, the kinetic energy showed significant sensitivity to smoothing parameters. Figure 13a shows the kinetic energy measured over time for each smoothing option and the smallest integration volume. All kinetic energy results shown are normalized ( KE * ) by the slug kinetic energy E = 1 8 D 2 o LU 2 p . Increasing smoothing caused drops of up to 30% of the slug kinetic energy, similar to what was observed for the Hill's vortex cases at lower resolutions (Fig. 9).
The kinetic energy increased with integration volume, as shown for the first four timesteps in Fig. 13b. Supplementary figure S7 shows each volume and smoothing option over all timesteps. While the total energy varies with time due to the inherent unsteadiness of the turbulent ring and the resolution of the experiment, the slopes with respect to integration volume size show the same trend at each time. The consistent increases in energy for volumes from 50 to 120 cm 3 suggest that a larger integration volume is required to completely capture the kinetic energy than was required to completely capture the impulse (Fig. 12a). For the spatiotemporally smoothed data, the energy levels off between the two largest integration volumes. However, for the Gaussian-smoothed data, the energy increases for the largest integration volume, suggesting a spike due to noise that was suppressed by the spatiotemporal smoothing. The energy increase from the 120 cm 3 volume to the 140 cm 3 volume also varies more between timesteps. Since changes in total energy as the integration volume expands should be consistent at neighboring times, it may be possible to detect energy increases from noise by looking at variations in profile shapes with respect to increasing integration volume.
To compare in-and out-of-plane resolutions, we compared the energy contributions of the X-and Z-velocity components. The strong axial jet of the vortex ring meant that significant kinetic energy in the volume was also due to Y-velocity components. The kinetic energy contribution of the Z-velocity components was greater than the energy contribution of the X-velocity components (Fig. 13c). On average, the Z-velocity components were 1.4× more energetic, greater than the 1.2× difference in interrogation window size. If we assume that the X-and Z-contributions should be equal, this is equivalent to a 6% ± 3% overestimate for the Gaussian-smoothed data and a 5% ± 1% overestimate for the spatiotemporally smoothed data, not including smoothing or domain selection effects. The kinetic energy of the Z-velocity components also fluctuated more between times than the X-velocity components.
Combined, these observations suggest interactions between the increased uncertainty and noise in the Z-direction, additional smoothing due to elongation of particles, and in this particular experiment larger interrogation windows in the Z-direction. These effects are likely difficult to decouple without prior knowledge of the flow field, highlighting a challenge for kinetic energy estimates in freely swimming scenarios where the orientation of a flow feature may change even between trials. These results also suggest that experiments be designed to avoid axial jet alignment with the out-of-plane (Z) direction to reduce uncertainty sources on kinetic energy, wake power, or efficiency estimates.

Impulse
Impulses for both window overlaps for the turning fish are shown in Fig. 14. We use a weighted centroid of vorticity as the impulse origin. Results were transformed to a coordinate system parallel and perpendicular to the final direction of fish travel ( −20 • rotation about the Y-axis). The increased window overlap has a smoothing effect on the time-resolved data, even though no temporal smoothing was applied. This smoothing is most noticeable when the propulsive vortex is shed and the fish body is mostly clear of the integration volume (t > 0.1 s). The impulse was noisiest at the beginning of the turn when the body was closest to the preparatory vortex, which was also located near the depth limits of the measurement volume.
In the direction of fish motion, there is a clear change and subsequent plateau in impulse; the initial and final impulses could be used for a stroke-averaged force estimate or differentiated with smoothing for instantaneous forces. Toward the end of the sequence, the impulse calculation captures contributions of subsequent tailbeats in the direction perpendicular to the fish's motion. All cumulative impulse fluxes were below 2 gcm/s (Supplementary Figure S9), suggesting that significant impulse losses due to the wake advecting outside the measurement volume did not occur over the timescale analyzed.
The values observed are comparable to previous 2D and 3D PIV studies on the same species of fish, which used circulation-based models to determine the impulse. Epps and Techet (2007) reports impulses of 32 gcm/s (for the preparatory vortex) to 112 gcm/s (for the propulsive vortex) for a rapidly turning giant danio. Mendelson and Techet (2015b) reports lower impulses of 7.7-11.9 gcm/s for forward swimming danio and 18.3-43.9 gcm/s for a slower turn.
Using a planar slice of the wake through the Y-position of the caudal peduncle (Supplementary Figure S8) and an axisymmetric vortex ring model (Eq. 3), we estimated impulse magnitudes of 68 gcm/s (50% overlap) and 59 gcm/s (75% overlap) for the propulsive vortex (calculation procedure in Supplementary Information). In comparison, the magnitude of the in-plane impulse change between t = 0.058 s and t= 0.138 s (using only parallel and perpendicular components since the planar estimate ignores I y ) was 57 gcm/s (50% overlap) and 67 gcm/s (75% overlap). Agreement could likely be further improved through tighter domain selection for the volumetric calculation or temporal filtering to reduce impulse fluctuations prior to t = 0.058 s.

Kinetic energy
Since no derivatives are taken in the calculation of kinetic energy, window size effects on kinetic energy measurement are entirely dependent on the smoothing filters used. Since smoothing was comparable between overlaps (kernel size and weights), the total kinetic energy over time was also similar (Fig. 15). As with the impulse, higher window overlap improved consistency over time, even without temporal filtering, and did not lead to lower energy values. Peak kinetic energies, capturing the fastest flow produced by the tail, were observed prior to the release of the propulsive vortex (t = 0.02-0.05 s). The higher fluctuation in energy over time in this interval also suggests increased noise due to the presence of the body. These factors highlight timing tradeoffs in kinetic energy analysis. The energy may be underestimated if evaluated after the body is completely clear of the wake, which also propagates into wake power and efficiency estimates. However, not all of the kinetic energy in the flow field at the peak time is necessarily lost to wake formation if the body is still interacting with the vortex. The wake energy was more constant once the propulsive vortex had been shed ( t > 0.1 s), but dissipation was a more significant consideration on this timescale. Dissipation rates calculated from the velocity gradients ranged from 0.014-0.034 mW for the 50% overlap case and 0.017-0.041 mW for the 75% overlap case. Cumulative dissipation was also slightly higher in the 75% window overlap case (0.0048 mJ versus 0.0041 mJ), roughly 10% of the peak kinetic energy in the volume. Dissipation was also likely still underestimated due to the limited resolution of nearbody boundary layer regions and is even more significant for lower Reynolds number organisms than the fish considered here. All energy losses due to fluxes out of the measurement volume (Eq. 10) were below 0.001 mJ (see Supplementary Figure S10), suggesting that while dissipation should not be neglected for this case, energy fluxes could be.

Conclusions and recommendations
Through these characterizations, we show that hydrodynamic impulse (Eq. 1) and wake kinetic energy (Eq. 5) can both be simply and accurately derived from volumetric velocity and vorticity field data. These results support increased quantitative use of volumetric velocimetry methods to assess swimming performance. Since the impulse and wake kinetic energy have extensive prior use in both 2D and 3D literature, these results also support increased comparisons between organisms and across studies.
Our results suggest that direct calculation of hydrodynamic impulse (Eq. 1) is particularly robust to noise and variations in processing and experiment design parameters.
The resolution needed to achieve an accurate impulse measurement is also achievable for volumetric PIV or PTV experiments. Furthermore, uncertainty in a direct impulse calculation from the vorticity field is comparable to uncertainty from determining the impulse using a vortex ring circulation model (Eq. 3), despite the additional origin selection variable, and may be more reliable for data with significant unsteadiness such as the turbulent vortex ring case.
The resolution needed to achieve an accurate kinetic energy measurement (Eq. 5) is also achievable for an idealized case, but this quantity is much more sensitive to setup, noise, and processing parameters. Energies calculated from volumetric velocimetry measurements are also particularly sensitive to trade-offs between noise and smoothing, which can arise especially in near-body regions of the flow field, as illustrated by the turning fish case. Kinetic energies should be calculated from velocity data with minimal post-processing if an upper bound on energy lost to wake formation is required (such as to evaluate efficiency). Any further postprocessing can cause underestimation, especially at limited spatial resolution. For cases such as the fish considered in this work, nontrivial dissipation rates suggest that wake kinetic energy is also a time-sensitive measurement. To compare energetics between different organisms or studies, compensation for differences in experimental parameters may be required for the most accurate comparisons.
These comparisons also suggest opportunities and outstanding questions in using volumetric velocimetry to evaluate swimming performance. In this work, we have tested simple smoothing methods to characterize post-processing effects; post-processing based on physical constraints such as the divergence-free filters used for fish wakes by Tu et al. (2022) may lead to increased accuracy. Approaches that combine metrics such as wake hydrodynamic impulse, which is both robust and widely reported in the previous literature, with capabilities for pressure estimates (as demonstrated for fish by Tu et al. (2022)) or full force and moment predictions (e.g., David et al. 2009;Rival and van Oudheusden 2017;Kang et al. 2017;Limacher et al. 2020) may also increase understanding of instantaneous force production during aquatic locomotion.