Calculating Energy Spectra from Drifters

The relations between the kinetic energy spectrum and the second order longitudinal structure function in two dimensions are derived, and several examples are considered. The forward conversion (from spectrum to structure function) is illustrated first with idealized power law spectra, representing turbulent inertial ranges. The forward conversion is also applied to the kinetic energy spectrum of Nastrom and Gage (1985), derived from aircraft data in the upper troposphere, and the result agrees well with the longitudinal structure function of Lindborg (1999), using the same data. The inverse conversion (from structure function to spectrum) is tested with data from 2D turbulence simulations. When applied to the theoretical structure function (derived from the forward conversion of the spectrum), the result closely resembles the original spectrum, except at the largest wavenumbers. However the inverse conversion is much less successful when applied to the structure function obtained from pairs of particles in the flow. This is because the inverse conversion favors large pair separations, which are typically noisy with particle data. Fitting the structure function to a polynomial improves the result, but not sufficiently to distinguish the correct inertial range dependencies. Furthermore the inversion of non-local spectra is largely unsuccessful. Thus it appears that focusing on structure functions with Lagrangian data is preferable to estimating spectra.


Introduction
Large scale turbulence in the ocean is responsible for the fluxes of energy, enstrophy and tracers across scales.Thus turbulence facilitates the dissipation of energy injected by the winds [1].However the exact nature of the turbulence remains uncertain because it is difficult to quantify kinetic energy spectra in the ocean.Satellite data, which has the spatial coverage required for such quantification, is limited in horizontal resolution, which is of order 100 km ( [2,3]).Shipboard current measurements can produce finer resolution estimates, down to several kilometers, in targeted regions ( [4,5]).But there are still relatively few studies in which spectra have been calculated, and the results often differ with region.
Relative dispersion concerns the separation of pairs of particles in flows.It was originally studied as a means of describing the spread of a passive tracer, beginning with Richardson's observations of smoke spreading from factory stacks [6].However it was later recognized that pair dispersion in a turbulent flow depends on the kinetic energy spectrum.Richardson's observations were later shown to be consistent with having a Kolmogorov turbulence spectrum [7].As such, pair dispersion offers a way of deducing kinetic energy spectra in geophysical systems.The advantage is that with particle pairs, one can in principle sample a broad range of scales.With the current generation of GPS-tracked surface drifters, one can resolve scales down to 10 m, two and four orders of magnitude smaller than is currently possible with shipboard measurements or satellite altimetry, respectively.
A number of studies have examined the relation between relative dispersion and turbulence spectra [8][9][10]."Local dispersion" occurs when pairs are influenced by eddies with scales comparable to the pair separation and "nonlocal dispersion" when they are dominated by larger eddies.Under local dispersion, the mean square pair separation (the dispersion) exhibits a power law growth, with an exponent related to the slope of the energy spectrum [11].Under non-local dispersion, the dispersion grows exponentially in time [8,11,12].
In the atmosphere, exponential dispersion was observed at separations below 2000 km in balloon experiments [13,14].Exponential dispersion has also been observed at the ocean surface, at scales below the Rossby radius of deformation [15][16][17].The balloon-and drifter-derived results agree with independent Eulerian estimates at these scales, in the atmosphere [18] and ocean [4].But pair dispersion estimates are often noisy and conflicting conclusions have been drawn.Some have suggested the aforementioned balloon dispersion is more indicative of local dispersion [19], while others have found evidence for local dispersion in the ocean, both at the surface [20,21] and subsurface [22,23] .
In all such cases though, the energy spectra were inferred from dispersion.However under certain assumptions the energy spectra can be calculated directly from velocity differences from pairs of drifters.The following focuses on this calculation in 2D flows, which are appropriate for large scale motion in the atmosphere and ocean.We examine the results in different cases, from idealized spectra to 2D turbulence simulations.

Spectra and structure functions
In a homogeneous flow, the two-point velocity autocorrelation and the kinetic energy spectrum are related via the Fourier transform [24]: where: R ij ≡< u i ( x + r)u j ( x) > , the brackets indicating an area average, over x.Similar relations apply in three dimensions as well, but we focus on two dimensions hereafter.For an isotropic flow, the spectrum depends only on the magnitude of the wavenumber, K, and the correlation is: Thus R ij (r) is the Hankel transform of Φ ij .The kinetic energy is just the trace of R ij , at zero separation: with summation implied for repeated indices.Using the orthogonality condition for Hankel transforms, the inverse relation is: In homogeneous and isotropic flows, the two point correlation can be determined from either Eulerian or Lagrangian measurements.
The velocity correlation in turn can be written in terms of the longitudinal ( f ) and transverse (g) velocity correlations [24]: where δ ij is one for equal indices and zero otherwise.Imposing continuity in 2D requires: [see also [25][26][27].Thus R ii can be written: Substituting this into (4) yields: after integration by parts and assuming that f vanishes at infinity.The inverse relation is: So E/(K 2 r) and K 2 r f are also a Hankel transform pair.Expressions ( 8) and ( 9) are the two dimensional counterparts of Batchelor's expressions 3.4.16[24].
A common measure used in analyzing turbulent flows is the second order velocity structure function: [ 24,28].In a homogeneous and isotropic flow, the longitudinal component, S 2l (r), is easily measured from the separation velocity of pairs of particles: This in turn is related to the longitudinal velocity correlation: where E l is the longitudinal kinetic energy, half the total energy in an isotropic flow.Thus: Expression ( 13) is similar to equation 3.1 of [11], except that the latter relates instead to the full structure function.What differs is the factor in parentheses (the full expression has (1 − J 0 (Kr)).The two expressions behave similarly however, acting as high pass filters (Fig. 1).In the limit of small separations, the filter in ( 13) is: This is half the limit for the full structure function, as the longitudinal component is half the total in an isotropic flow.For larger separations the function approaches one, and S 2l (r) asymptotes to the twice total kinetic energy.The energy can be determined from the inverse relation, which follows from (8): This also involves a high pass filter, K 2 r 2 J 1 (Kr).However this has no asymptotic limit for large separations; this complicates the inversion, as seen hereafter.Nevertheless, a particle-derived structure function can in principle be used to calculate the energy spectrum.

Idealized spectra
We'll examine the forward conversion first, calculating the structure function from the spectrum.Consider an idealized spectrum: The spectrum has a finite range (rather than extending from K = 0 to K = ∞) and thereby avoids associated singularities.Using the asymptotic limits for the high pass filter [11], the structure function ( 13) is approximately: for separations such that K 0 < 2π/r < K 1 .If α > 3 (a non-local spectrum) and if the lower wavenumber, K 0 , is small enough, the integral is approximately: The structure function increases as r 2 , regardless of how steep the spectrum is.If the slope is such that 1 < α < 3 (local spectrum) and K 1 is sufficiently large, the second term dominates: For a Kolmogorov spectrum (α = 5/3), the structure function scales as r 2/3 .However unless the spectral range is long (K 1 K 0 ), the other terms will contribute and cause deviations.For separations such that r > 2π/K 0 , the structure function is constant: Thus the lower wavenumber limit, K 0 , determines the scale at which S 2l flattens out.
The result with α = 3 is shown in Figure (2).The case discussed above is plotted in blue in the left panel, with K 0 = 1 and K 1 = 256.The two other spectra peak instead at K = 5.The structure functions in the right panel were calculated numerically from (13).
The structure function for the spectrum peaked at K 0 is shown in blue in the right panel.This increases somewhat slower than r 2 , illustrating the effect of having a finite wavenumber range.The two other spectra exhibit a similar increase at small scales but level off sooner, due to having a larger peak wavenumber.Note too that the difference at the low wavenumbers between the two spectra with K = 5 has little impact on the structure functions, due to the high pass filter in (13).
One obtains nearly identical results if one replaces the high pass filter in (13) with the function: So the transform is relatively insensitive to the form of the filter.Furthermore, the results with local spectra are qualitatively similar to those seen above.

The Nastrom and Gage (1985) spectrum
For a more realistic example, we consider energy spectra from the upper troposphere.Nastrom and Gage calculated zonal and meridional wavenumber spectra from data collected from commercial aircraft over the United States [18].Lindborg calculated longitudinal and transverse structure functions, using the same data [26].Lindborg rationalized the structure functions using asymptotic limits appropriate for the two inertial ranges seen in the Nastrom and Gage spectra, specifically K −3 for separations from roughly 2000 km to several hundred kilometers and K −5/3 at smaller scales.However, one can estimate the longitudinal structure function directly from the spectra, using (13).We extracted the zonal and meridional spectra from Figure 3 of [18], 1 and then used the total energy to evaluate (13), numerically.
The resulting second order structure function is shown in the right panel of Figure (3), with the longitudinal structure function of [26] indicated by the dotted curve.The dashed line indicates twice the total energy.The two curves agree well, except at the smallest scales where the predicted S 2l is less than the observed.This may reflect deviations from isotropy and/or 2D non-divergence at small scales [e.g.29].Nevertheless the overall agreement is good.

Turbulence spectra
We now perform the inverse relation, calculating the spectrum from the structure function.For this we employ data from several 2D turbulence simulations, studied previously by [17].The numerical code [30][31][32] solves the barotropic vorticity equation: where ψ is the velocity streamfunction, such that: is the Jacobian operator and F and D are the applied forcing and dissipation.The domain is doubly-periodic, with 512 2 grid points.The forcing is isotropic and applied in specified wavenumber ranges, with random phases in space and time.The amplitude was set so that the final (dimensionless) kinetic energy was approximately 1.0.For dissipation, we employed linear (Rayleigh) friction, The model also has an exponential cut-off filter which removes energy at the smallest scales [33].
A snapshot of the streamfunction in the equilibrated state of one run is shown in the left panel of Figure (4).The forcing in this case is applied in the wavenumber range K = [30,35] and the Rayleigh coefficient was R = 0.1.There are numerous large eddies, but smaller scale features are seen as well.The large eddies exceed the forcing scale (π/30 ≈ 0.1), suggesting an inverse energy cascade [e.g.34].
When the stationary state was reached, 2000 particles were deployed on a uniform grid and advected with a fourth-order interpolation scheme.The initial particle separation was r 0 = 0.01, comparable to the model grid size.The trajectories are shown in the right panel of Figure (4).
The time-average energy spectrum was shown in Figure 9 of [17] and is reproduced here in the left panel of Figure (5).The spectrum is steep from the forcing range (K = [30,35]) to roughly K = 200, where the exponential filter is acting.The slope is nearly −4 and is thus greater than the theoretical value (−3) for a 2D enstrophy range [35].This is due to the Rayleigh damping; simulations without this produce a slope closer to −3, as seen below.The spectrum is nevertheless non-local above the forcing wavenumber, implying the same impact on pair separations as in an enstrophy range [11].There is also a short energy range, with a slope of K −5/3 , from the forcing wavenumbers up to K = 10.This is in line with an inverse energy cascade, as inferred from Figure (4).
The second order longitudinal structure functions are shown in the right panel of Figure (5).Three curves are plotted.One (blue solid) is that predicted from the energy spectrum using (13).The second (red dashed) was obtained using velocity differences from all particle pairs in the flow.The Peer-reviewed version available at Fluids 2016, 1, 33; doi:10.3390/fluids1040033third (yellow dash-dot) was calculated using velocities only from pairs deployed together, "original pairs" in the terminology of [13].
The three curves are very similar.All increase nearly as r 2 from the grid scale (r = 0.01) up to the forcing scale (r = 0.1).The curves flatten out above r = 1.0.The r 2/3 dependence expected for an energy range is not seen, due most likely to its limited wavenumber extent.Furthermore, the two particle-derived estimates are very similar.This implies that the choice of original or "chance" pairs does not yield qualitative differences in S 2l ; it merely effects the number of samples (there are many more chance pairs).
To test the inverse operation, we first use the structure function calculated from the energy spectrum.We do this numerically, using (14). 2 This is of course just a test of the invertibility of the Hankel transforms linking the two measures (since the theoretical structure function derives from the original spectrum).The result is shown in Figure (6).The reconstruction is good in the wavenumber range below the forcing range, but at higher wavenumbers the predicted spectrum deviates from the original substantially, even becoming negative at some wavenumbers.
The failure at large wavenumbers is due to the filter function in ( 14), K 2 r 2 J 1 (Kr).As seen in Figure (7), this oscillates and grows with separation.Thus the small wiggles present in the structure function where it levels off (Fig. 5) affect the energy spectrum at large wavenumbers.
However, an alternate approach is possible.The longitudinal velocity correlation, f (r), is easily obtained from the second order structure function, as in (12).Then one can calculate the two point correlation, R ii (r), using (7), and the energy spectrum then follows from (4).This is preferable to (14) because the filter in ( 4) is KrJ 0 (Kr).While this also increases at large separations, it does so more slowly than K 2 r 2 J 1 (Kr) (Fig. 7).
The spectrum obtained using this alternate approach is shown in Figure (8).Now a greater range of wavenumbers is reproduced, up to roughly K = 100.Hence this approach produces a more satisfactory inversion, and we use it hereafter.
Next we turn to the particle-derived structure functions.For this, we use the S 2l obtained from all particle pairs (the results with the original pairs are similar but noisier).The reconstructed spectrum is very noisy, as indicated by the red dots in Figure (9).This again is due to the inversion (4) favoring contributions from large separations.
The results are improved somewhat if one fits the the structure function to a polynomial, to remove the oscillations at large separations.We did this by fitting the logarithm of S 2l with a sixth order polynomial, as shown in the left panel of Figure (10).As seen, the behavior at small scales is well captured by the fit while the wiggles at large scales are greatly reduced.
Shown in the right panel are the two point correlations, R ii (r), from the particle-derived S 2l (in red) and from the polynomial fit S 2l (in yellow).Also plotted is the two point correlation obtained from the theoretical S 2l (in blue).The latter decreases from a maximum at r = 0.01, crosses zero near r = 0.9 and then increases gradually towards zero. 3The particle-derived correlation is similar at small separations but is much noisier at larger separations, oscillating between positive and negative values.The polynomial fit remedies this somewhat by removing the oscillations.However, the resulting correlation still differs from the theoretical prediction.
Using the two point correlation from the polynomial fit in (4) yields the yellow dotted curve in Figure (9).The spectrum is still noisy, particularly at higher wavenumbers, but is significantly improved over the estimate using the original structure function.That said, one is not able to detect the two inertial ranges present in the original (Fig. 5).The calculated spectrum instead rolls off at higher wavenumbers.
2 [36] discuss the numerical evaluation of Hankel transforms.Here we use transforms which are available in Matlab.

3
The correlation is adjusted so that the last value is zero, as expected for the Hankel transform [36].The results from two other simulations are shown in Figures (11)(12).The first is for a pure energy cascade, forced in the range K = [100, 120], with a Rayleigh coefficient of R = 0.1.Plotted are the original spectrum (in black), the reconstruction from the theoretical S2 (in red), from the particle S2 (green dots) and from a 6th order polynomial fit of the particle S2 (blue).The spectrum exhibits a well-defined K −5/3 inertial range from the forcing range down to K = 10.The spectrum from the theoretical S2 matches the original to roughly K = 100, as in the previous case.The spectrum from the raw S 2l is again noisy and inconclusive, while that from the polynomial fit is closer to the actual spectrum.However the latter is also steeper than the true spectrum, giving an incorrect estimate of the slope in the inertial range.
The results with a pure enstrophy cascade, forced in the range K = [1, 5] and with no Rayleigh damping (R = 0), are significantly worse (Fig. 12).The slope in the inertial range is close to, though slightly steeper than, K −3 .However the spectrum from the theoretical S2 (red dotted curve) is accurate only at the smallest wavenumbers.It begins to oscillate around K = 10 and these oscillations increase with wavenumber; the dotted curve indicates the upper branch of the curve, and the lower values are actually negative.The spectrum derived from the polynomial fit (blue dotted curve) deviates from the original spectrum at even smaller wavenumbers, and exhibits negative values even below K = 10.The upper branch mirrors that from the theoretical S2, but the values are roughly 10 times larger.The spectrum from the raw data (green dots) is again poor.
As the curves are oscillating between positive and negative values, we tried averaging the values in logarithmically-spaced bins.This produced estimates closer to the observed spectrum (not shown).But the results remain fairly poor.
Thus the inverse conversion essentially fails for non-local spectra.Other methods of fitting the structure function were attempted, and the two point correlation was also fit.Generally the polynomial fit described above yielded the best results.It may be that more refined methods of filtering the structure functions will be more successful, but the differences between the particle-derived correlation and the theoretical one (Fig. 10) suggests this is unlikely.

Summary and discussion
Relations between the kinetic energy spectrum and the second order velocity structure function were derived in two dimensions.These are related to those given by [11] but are specifically applicable to longitudinal structure functions, which are easily measured with pairs of particles.The relations involve Hankel transforms, and can be evaluated with routines found in Matlab.
Structure functions derived from idealized spectra with power law dependencies exhibit power law dependencies on separation, as described previously [e.g.8,11].However the dependencies obtained with inertial ranges of finite extent differ somewhat from familiar relations which assume infinite ranges.This probably explains why the r 2 behavior is rarely observed in situations where the spectra are believed to be non-local [e.g.37].We also applied the transform to the kinetic energy spectrum in the upper troposphere of [18].The resulting second order structure function closely resembles the longitudinal structure function calculated from the same data by [26].
The relations were then applied to data from several 2D turbulence simulations.First the structure function was calculated from the spectrum using (13), then the inverse relation (14) was used to recover the energy spectrum.The spectrum was reproduced at small wavenumbers but not at larger wavenumbers, because (14) has a filter which magnifies contributions from large separations, hindering the Hankel transform.However, calculating the two point velocity correlation first and then using (4) produced a spectrum which matched the original over most of the wavenumber range.
The problems associated with large separations are amplified when using particle-derived structure functions, which are noisy precisely at the largest scales.The spectral estimates are improved if one fits the structure functions with a high order (6th in this case) polynomial, but they deviate nonetheless from the original spectra.The results also are better when the spectra are local (shallower than K −3 ).The transform in its present form functions poorly with non-local spectra.
However it may be that application in future of more sophisticated data processing will improve this situation.
Thus the forward transform (spectrum to structure function) is more robust than the inverse relation.This follows from the difference in the high pass filtering inherent in the two operations.While the filter asymptotes to one for large separations with the forward transform (13), it increases with separation for the inverse (4).A similar issue pertains to the relations in three dimensions, as seen in equations 3.4.16 of [24].As such it is easier working with the forward transform to assess a given structure function, like that of Lindborg's, than to do the inverse.
The present results imply that it is probably better to examine structure functions from Lagrangian data then to extract the corresponding energy spectra.As seen here, the structure functions are fairly robust, agreeing well regardless of whether chance or original pairs are used.Nevertheless, structure functions integrate contributions across wavenumbers and thus are not as "clean" as energy spectra.
Several recent studies have examined decomposing energy spectra into rotational and divergent components [29,38,39].This usefully separates contributions from large scale balanced motions (e.g.vortices) and smaller scale phenomena like waves.A similar decomposition can be applied to structure functions [40,41].The present relations can be used to test consistency between the two, as shown here in relation to the atmospheric spectra.

Figure 1 .
Figure 1.The high pass filter, used in converting the energy spectra to a second order structure function.The filter of Bennett (1984) is shown for comparison.

Figure 2 .
Figure 1.The high pass filter, used in converting the energy spectra to a second order structure function.The filter of Bennett (1984) is shown for comparison.

Figure 3 .
Figure 3.The zonal velocity spectrum from [18] (left panel).The corresponding 2D longitudinal structure function is shown by the solid curve in the right panel.Also shown is the structure function (dotted curve) calculated from data, from [26].The latter has been divided by 2. The dashed horizontal line is the twice the zonal energy, calculated from integrating the zonal spectrum.

Figure 4 .
Figure 4. Streamfunction (left panel) and particle trajectories (right panel) from a two dimensional turbulence simulation forced in the range K = [30, 35].

Figure 5 .Figure 6 .
Figure 5. Kinetic energy spectrum (left panel) from the turbulence simulation, reproduced from Fig.9of[17].The structure functions (right panel) calculated from the energy spectrum (solid curve), from all particle pairs (dashed red) and from pairs deployed together (dashed yellow).

Figure 7 .Figure 8 .
Figure 7.The "filter" function used in the spectrum reconstruction in (14) (dashed curve).Also shown is the function used in the alternate reconstruction in (4) (solid curve).

Figure 9 .Figure 10 .
Figure 9. Kinetic energy spectra using the reconstruction based on (4), for the particle-derived second order structure function (red dots).Also shown is the reconstruction based on a 6th order polynomial fit of the structure function (yellow dotted curve).

Figure 11 .Figure 12 .
Figure 11.Spectral reconstructions from a simulation of a 2D inverse cascade, forced in the range K = [100, 120].The original spectrum is in black and the reconstruction from the theoretical structure function is in red.The green dots show the spectrum from the raw particle S 2l and the blue curve shows that obtained with the polynomial fit S 2l .