Second order structure functions for higher powers of turbulent velocity

We experimentally study the temporal second-order structure functions for integer powers of turbulent fluid velocity fluctuations , in three dimensional (3D) and two dimensional (2D) turbulence. Here is a composite time-series constructed by averaging the concurrent time-series () sampled at N spatially distributed Eulerian points. The N  =  1 case has been extensively studied for velocity fluctuations (m  =  1) and to a lesser extent for m  >  1. The averaging method in case of N  >  1 diverges from the Kolmogorov framework and has not been studied because fluctuations in are expected to smooth with increasing N leaving behind uninteresting large-scale mean flow information, but we find this is not so. We report the evolution of scaling exponents for in going from a single (N  =  1) to a spatial average over several Eulerian points (). Our 3D experiments in a tank with rotating jets at the floor show for all m-values in agreement with prior results and evolves to an asymptotic value of . The evolution of follows the functional form , where points is the only fit parameter representing the convergence rate constant. Results for the 2D experiments conducted in a gravity assisted soap film in the enstrophy cascade regime are in sharp contrast with their 3D counterparts. Firstly varies polynomially with m and asymptotes to a constant value at m  =  5. Secondly, the evolution of is logarithmic , where A and B are fit parameters and eventually deviates at large N and asymptotes to for all m. The starkly different convergence forms (exponential in 3D versus logarithmic in 2D) may be interpreted as a signature of inter-scale couplings in the respective turbulent flows by decomposing the two-point correlator for into a self-correlation and cross-correlation term. In addition to aiding in the theoretical development, the results may also have implications for determination of resolution in 2D turbulence experiments and simulations, wind energy and atmospheric boundary layer turbulence.


Introduction
The Kolmogorov theory of 1941 (K41) [1][2][3] holds a special place among all statistical hydrodynamic approaches [4,5] to turbulence. The first of Kolmogorov's two works from [1] (K41A) concerns the second moment of longitudinal velocity differences (∆u || (r) ≡ u || ( R + r) − u || ( R)) or the secondorder structure function, which predicts (∆u || (r)) 2 ∼ (εr) 2/3 . Here ε is the mean energy flux per unit mass and η < r < l 0 denotes a length scale within the inertial range bounded between the dissipative (η) and integral (l 0 ) scales. Whereas a proper experimental deduction of this result that is faithful to theoretical requirements demands an instantaneous field We experimentally study the temporal second-order structure functions for integer powers of turbulent fluid velocity fluctuations S m 2 (N, τ ) ≡ (U m (t + τ ) − U m (t)) 2 , in three dimensional (3D) and two dimensional (2D) turbulence. Here U m (t) = (1/N)Σ N i=1 u m i (t) is a composite time-series constructed by averaging the concurrent time-series (u m i (t)) sampled at N spatially distributed Eulerian points. The N = 1 case has been extensively studied for velocity fluctuations (m = 1) and to a lesser extent for m > 1. The averaging method in case of N > 1 diverges from the Kolmogorov framework and has not been studied because fluctuations in U m are expected to smooth with increasing N leaving behind uninteresting large-scale mean flow information, but we find this is not so. We report the evolution of scaling exponents where A and B are fit parameters and eventually deviates at large N and asymptotes to ζ m (N 1) = 2.0 ± 0.1 for all m. The starkly different convergence forms (exponential in 3D versus logarithmic in 2D) may be interpreted as a signature of inter-scale couplings in the respective turbulent flows by decomposing the two-point correlator for U m into a selfcorrelation and cross-correlation term. In addition to aiding in the theoretical development, the results may also have implications for determination of resolution in 2D turbulence experiments and simulations, wind energy and atmospheric boundary layer turbulence.
Keywords: turbulence, spectra, fluctuations (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. measurement, early experiments relied on Eulerian timedomain measurements [6], which then exploited Taylor's Hypothesis [7] to interpret data [8][9][10].
Understanding the spatio-temporal character of turbulence fluctuations however requires both space and time spectra [11], which albeit possible in simulations [12] is experimentally still difficult at least in three dimensional (3D) flows. Whereas spatial spectra from instantaneous snapshots do recover most of the information for statistically stationary flows and chart out the K41A predictions, the same is not necessarily true for Eulerian time-domain measurements [13]. Be that as it may, a temporal Eulerian probe still represents a justifiable physical measurement, as is indeed the norm in many an application from meteorological monitoring and atmospheric turbulence [14] to wind energy [15,16]; such signals contain information worthy of deduction in their own right. Sometimes the measured quantities are functions of higher powers of turbulent velocity, as in case of wind power [16], and involve higherorder turbulence spectra [17].
We report experiments in both three dimensional (3D) and two dimensional (2D) turbulence (2D is in enstrophy regime) with focus on temporal, Eulerian, second-order structure functions S m 2 (N, τ ) ≡ (U m (t + τ ) − U m (t)) 2 for higher integer powers of a composite velocity time-series ( U m ; m 1). The composite time-series U m (t) is constructed from an average over the concurrent, individual time-series (u m i (t)) sampled at spatially discrete, Eulerian points in the turbulent field and is defined as Here u m i (t) denotes the time-series for the mth power of velocity at the ith location.
Since an average over the N Eulerian points denotes a spatial average, U m (t) represents a composite time-series signal that is spatially averaged at each instant. Standard turbulence analysis follows averaging proposed in K41 and not the kind we describe for a simple reason. When an increasing number (N) of Eulerian points participate in the spatial averaging prior to computing the spectra, local fluctuations at smaller scales are expected to smooth out leaving behind only largescale mean flow characteristics in U m (t), representative of a coarse grained version of a field average in the N → ∞ limit. Accordingly, the quantity S m 2 (N 1, τ ) which captures information about fluctuations and their correlations is not expected to show any interesting scaling behavior, if at all it does exhibit any scaling, i.e. if S m 2 (N, τ ) = A m (N)τ ζm(N) , it is expected to result in A m (N) = 0 and ζ m (N) = 0, provided N is sufficiently large. In contrast, not only do our 2D and 3D experiments reveal scaling for ζ m (N), it evolves with N and asymptotes to a constant value. The evolution of ζ m (N) and its asymptotic convergence are markedly different between the 2D enstrophy sub-range and 3D. Firstly the convergence is exponential in 3D whereas it is logarithmic in 2D and secondly, the 3D asymptotic value converges to an m-dependent value of ζ m (N) = 2m/3 whereas in 2D it asymptotes to a constant value of ζ m (N 1) = 2.0 ± 0.1. The common starting point for our analysis and K41 framework is the case of velocity fluctuations (m = 1) at a single Eulerian point (N = 1); this is the most well studied situation where Taylor's [7] or random sweeping hypothesis [13,18] usually apply. Spectra for higher powers of velocity (m > 1) at a single Eulerian point (N = 1) too have been studied [17,19] and are referred to as 'Higher-order Spectra' since their analysis was in Fourier domain. Our analysis being in the timedomain, we use the term 'spectra' interchangeably to denote temporal second-order structure functions owing to their equivalence and to avoid cumbersome use of 'second-order structure functions for higher powers of turbulent velocity' at each instance. Higher-order spectra should not be confused with higher order structure functions of velocity fluctuations, the two quantities are quite distinct. To clarify further on this, for simplicity, let us consider velocity fluctuations (u 1 ) at a single Eulerian point (N = 1). Firstly, one can compute the higher order structure functions of arbitrary order n defined as S 1 n (τ ) ≡ (u 1 (t + τ ) − u 1 (t)) n , but only the n = 2 case has a Fourier domain counterpart, namely the velocity spectrum. Structure functions of order n = 2, S 1 n =2 (τ ) have no equivalent spectral counterparts in the Fourier domain. Likewise, the higher-order structure functions can be generalized for arbitrary powers of velocity (u m ) giving the most general for mulation S m n (τ ) ≡ (u m (t + τ ) − u m (t)) n . We only study the second order structure functions of this general formulation for arbitrary powers of velocity, viz. S m 2 (τ ) ≡ (u m (t + τ ) − u m (t)) 2 , whose spectral counterparts in Fourier domain are the Higherorder spectra. Secondly, as stated above, it is owing to this equivalence between the Fourier domain quantities (Higher order spectra) and physical domain quantities (second-order structure functions for mth power of velocity) that we use the term 'spectra' interchangeably to avoid cumbersome language.
As soon as we move beyond N = 1 case and employ N > 1 points to construct U m (t), our method of analysis no longer conforms to K41 framework. Careful scrutiny of multi-scaling behavior in turbulence has established that scaling exponents are ever so sensitive to the type of, and the order in which statistical averaging is conducted when analysing turbulence data (see [20] and references therein). Then it remains to properly interpret the scalings and draw appropriate inferences to aid in theoretical development of our proposed statistical averaging scheme. Before delving further, it behooves us to explain the motivation to attempt this statistical analysis in the first place. Our analysis was motivated in part by recent work on wind power fluctuations [16]. In the wind power problem, the power P(t) generated by a turbine varies with wind speed as P(t) ∼ u 3 (t). Explanation of the turbine-scale wind power fluctuation spectrum naturally takes one into a nuanced analysis of higher-order turbulence spectra, i.e. N = 1, m > 1 case. Furthermore, the aggregate wind power feeding the electrical grid is derived from turbines located in geographically distributed wind plants. Understanding the fluctuations in aggregate wind power at the grid level therefore requires summing individual time-series of spatially distributed Eulerian points (turbines) into a composite time-series signal, and only then computing the aggregate wind power fluctuation spectrum. The questions we sought to address in this experimental study were three fold. Firstly, what is the asymptotic field average value of ζ m (N) as a function of m? Secondly, what is the convergence rate of the exponent ζ m (N) to the asymptotic field limit as a function of N? Finally, what is the dependence of ζ m (N) and its convergence to field limit on flow dimensionality, since it is known that 2D turbulence has it's distinct complexity through presence of two cascades?
The rest of this article is structured as follows: sections 2 and 3 detail the experiments for 3D and 2D turbulence, respectively. The results are discussed in section 4 and we conclude with a summary in section 5.

Background
K41A [1] relates the second moment of longitudinal velocity differences, an experimentally measured quantity, to ε through a relation that scales self-similarly with an inertial range spatial separation r: Here (1) carries the physical interpretation of small eddies nesting within larger eddies in a self-similar structure within the inertial range of three-dimensional, isotropic, incompressible, homogeneous turbulence. Henceforth, we drop the subscript || denoting the longitudinal velocity component for sake of brevity and use u ≡ u || instead.
For an experimental determination of equation (1) to be physically consistent with K41A, one must simultaneously measure longitudinal velocity differences at various spatial separations r ≡ |( R + r) − R| at several spatial locations within the turbulent field. In other words, both u( R + r) and u( R) must be measured at the same time instant. Owing to experimental limitations, early turbulence experiments involved Eulerian (spatially fixed) point measurements over a time period t. Then invoking Taylor's hypothesis (TH) or the 'frozen turbulence assumption' [7] permits a linear transformation r ≡ uτ , where τ is an inertial range time scale corresponding to the length scale r. This in turn recasts equation (1) in time domain as: TH applies when the mean flow speed (u) is large relative to fluctuation magnitude (RMS velocity intensity I = u rms /u << 1 so that the eddies transported past the probe do not measurably evolve over the measurement duration. An eddy of size r in the inertial range has a fluctuation time scale τ = r/ (∆u(r)) 2 associated with it. When swept past an Eulerian probe by u, the eddy experiences a doppler shift and registers a time scale shorter than τ . This doppler shift is however uniform across eddies of all sizes since u (zero frequency mode) has no associated time scale of it's own. As opposed to TH, the random sweeping hypothesis (RSH) applies when there is no mean flow or if u rms ∼ u. This scenario was carefully considered by Kraichnan [13], Lumley [21], and Tennekes [18] among others [22,23]. Unlike TH, the RSH regime must contend with measurable eddy distortion over measurement duration, and u rms being non-negligble contributes to eddy transport past the Eulerian probe. As explained by Tennekes [18] in a physically intuitive manner, the integral scale eddy (size l 0 and velocity fluctuation magnitude u rms ) convecting a smaller eddy of size r nesting within it induces an oscillation in that eddy of time scale r/u rms and doppler shifts the eddy's fluctuation time scale τ to shorter time scale. This doppler shift is however scale dependent; smaller eddies with shorter τ experience a greater doppler shift than larger eddies. Consequently, RSH predicts the second-order structure function scales as: The replacement of u in equation (2) with u rms in equation (3) may seem straightforward if one naively interprets it as a case of u rms appearing in the spectrum because it convects small eddies in place of u. Yet, it leads to a physical situation fundamentally at odds with K41A assumptions. A dimensional analysis subject to constraints imposed by K41A assures us that ε is the only global parameter that links all spatial scales in the inertial range of turbulence. In particular, the small and large scales are statistically independent of each other and the small scales are not influenced by the presence of the large scales that convect them. Taylor's hypothesis does not contradict this picture because the mean velocity prefactor (u 2/3 ) in equation (2) is not linked to any inertial range length or time scale. However, the introduction of integral scale velocity fluctuation u 2/3 rms in equation (3) violates this assumption. All eddies are influenced by the integral scale eddy when it convects them away, and they experience a scale dependent doppler shift, or doppler broadening in RSH.
Unfortunately, it is very difficult to experimentally distinguish TH from RSH except for an ad hoc rule that u rms /u 1. Firstly both equations (2) and (3) predict the same τ 2/3 scaling. Secondly, the pre-factors (uε) 2/3 in equation (2) and (u rms ε) 2/3 in equation (3) are of same magnitude. This interpretational difficulty has historically shared close relation with developments in higher-order spectra, which we now briefly review.
To the best of our knowledge, Dutton and Deaven [19] were the first to consider the scaling of higher order spectra for algebraic (m) powers of turbulent velocity u m (see also discussion in [17]). Dutton and Deaven's calculation was performed in Fourier domain and essentially constituted a dimensional analysis, which we present below in the time domain. If K41A type inertial range scaling exists for second-order structure functions of integer powers of velocity (u m ,m > 1), i.e. S m 2 (r) ≡ (∆u m (r)) 2 are functions of ε and r alone, dimensional analysis shows: The above dimensional argument immediately yields α = β = 2m 3 . S m 2 (r) scales the same as S 1 2m (r) as a consequence of dimensional constraints arising from requiring S m 2 (r) be a function of only ε and r. In other words, under K41A requirements the quantity S 1 2m (r) ≡ (∆u 1 (r)) 2m scales as S 1 2m (r) ∼ (εr) 2m/3 but these are higher-order structure functions (of order 2m) for velocity (u 1 ) and are unrelated to S m 2 (r) ≡ (∆u m (r)) 2 , which too is expected to scale as S m 2 (r) ∼ (εr) 2m/3 due to dimensional constraints. Applying TH to equation (4) we obtain: Here C m are constants of undetermined order. Dutton and Deaven's measurements [19] however did not confirm the scaling expectation in equation (5). Instead they found S m 2 (τ ) ∼ τ 2/3 , m 1. More extensive measurements by Van Atta and Wyngaard [17] confirmed this for up to m = 9, which they complemented with theory in Fourier domain. We reproduce their arguments in time domain below. Consider the case of m = 2: Performing a simple decomposition of the statistical average gives to leading order: where, in equation (6c), we used the fact that S 1 2 (τ ) = u 2 rms − u(t)u(t + τ ) . The above result can be extended to arbitrary u m by utilizing the algebraic decomposition (a m − b m ) = (a − b)((m − 1) order term), which immediately gives: Simple substitution of equation (2) for TH or equation (3) for RSH in equation (7), then provides the scaling form for all m: We draw the reader's attention to the presence of u rms in equation (8) for both TH and RSH regimes. This large scale influence automatically enters higher-order spectra (m > 1) owing to the fact that whereas velocity differences are Galiean invariant, differences in higher powers of velocity do not satisfy Galilean invariance.
Thus far, we have considered spatial spectra for a turbulent field or time spectra at a single Eulerian point (N = 1) for m 1. Building upon the reviewed literature, we now consider several spatially distributed Eulerian points (N > 1), each of which registers a concurrent signal for an arbitrary power of velocity u m (t). Now if one were to sum these concurrent time-series from spatially distributed individual Eulerian points, whether their respective signals are correlated or not, into one composite signal, one obtains an averaging of the signal of the form: where each i-value represents an Eulerian point in the turbulent field. We are interested in learning what is the spectrum of the composite signal U m (t), i.e.: (10) Probing equation (10) as a function of increasing N allows us to track the evolution of S m 2 (N, τ ) ∼ τ ζm(N) towards the field averaged limit N → ∞, but in the time-domain thus setting it apart from the core K41A approach. Since no theoretical expectation is available for this deviation from K41A, we now proceed to describe the 3D experiments.

Experiment
The 3D experiments were conducted in a setup (figure 1) previously used to study free surface flows [24][25][26][27] and closely followed the measurement method described in [28]. The square tank of side length 1 m was filled with water to a height of 0.3 m. Turbulence was generated by means of 8 horse power pump, which removed water from side ports placed at the tank bottom and recirculated back through a system of circulating jets placed in a square grid on the tank floor. The turbulence generated by this injection scheme is not homogeneous in the vertical direction as it loses intensity with height, but the turbulence is homogeneous at any given horizontal plane within the bulk flow. This injection scheme also does not result in a mean flow, hence our 3D experiments fall squarely within the random sweeping regime; for this reason the turbulent intensity I is undefined.
A 5.5 W laser beam of 532 nm wavelength (coherent verdi diode pumped solid state laser) was passed through a cylindrical lens to generate a laser sheet and vertically positioned 4 cm below the free surface of water in the tank to illuminate a two-dimensional plane within the bulk fluid. Neutrally buoyant hollow glass spheres (TSI Inc. Model No. 10089) of nominal mean diameter 8-12 µm and density 1.05-1.15 g/cc were suspended in the fluid as tracer particles. A high speed camera (Vision Research Phantom v640) positioned above the tank recorded flow images of the laser illuminated plane over a square area of side length 5 cm at the center of the tank far from the walls. The camera collected images at 500 frames per second (fps) at a resolution of 1024 × 1024 pixels. The spatial and temporal resolution in our experiments are therefore 5 cm/1024 pixels = 48 µm/pixel and 1/500 fps = 2 ms/frame respectively. Particle imaging velocimetry (PIV) codes written in house were then employed to construct the velocity field from the raw image data. The tracer particle concentration limited the PIV grid resolution to 8 pixels. Each data set spanned a little less than 20 large eddy turnover times and 20 independent data sets were collected in quick succession under identical experimental conditions. The relevant parameters measured in the experiments are tabulated in table 1.
The computation of structure functions S m 2 (N, τ ), N 1 began with the PIV velocity vector field. Despite lack of mean flow, there was small anisotropy in the flow due to location of jets underneath the flow. We performed spherical harmonic decomposition to recover the isotropic sector from the PIV fields [29]. Surrogate fields, one for each m-value were generated by raising the velocity value at each grid point to the mth power. With increasing m-value, the fluctuations amplify quickly and overwhelm the statistics of u m , specifically the tails of the distributions. Our experimental data was sufficient for computations up to m = 4 in 3D.
Since each PIV field represents a snapshot for a given time step, selecting a single Eulerian grid point and tracking its velocity fluctuations in time across PIV snapshots provided a time-series for the N = 1 case. Since the PIV fields spanned a square window of side length greater than the correlation length, two random points were selected in each field ensuring their separation exceeded the correlation distance computed from the PIV field. From the recorded time-series spanning several correlation times across all data sets, the second-order structure functions for different m-values were retrieved by taking the difference ∆u m (τ ) = u m (t + τ ) − u m (t) for various τ -duration windows with a maximum window of τ 0 = l 0 /u rms , the turnover time as calculated. The S m 2 (N = 1, τ ) were then computed by squaring the difference (∆u m (τ )) 2 and averaging over several such differences across the time-series and all experimental runs to yield S m 2 (N = 1, τ ) = (∆u m (τ )) 2 . For N > 1, the same process was repeated, but each of the N points was selected randomly across the PIV grid, representing random inter-point distances, which could lie within or beyond correlation length. Then the additional step of spatial averaging of the individual time-series was performed to generate a single composite time series Using this composite time series U m (t), we then performed the procedure detailed above to compute S m 2 (N > 1, τ ) ≡ (U m (t + τ ) − U m (t)) 2 . This procedure was repeated for several N-values until all the grid points were covered. For each value of N, several randomly generated configurations of discrete Eulerian points were sampled, one of those configurations being the standard square PIV grid configuration. The results we quote are the exponents obtained from averaging over power-law scalings of S m (N, τ ) over all randomly generated configurations of N.
Once the structure functions were computed as per the above procedure, their scalings were analyzed for each value of m, and for each value of N. The scalings for S m 2 (N, τ ) ∼ τ ζm(N) were then obtained from a power-law fit to data to recover the scaling exponents ζ m (N). We note that the exponents ζ m (N) we present for 3D (figure 4) and 2D (figure 8) are average values, i.e. for each value of N, several randomly generated configurations were employed, each resulting in a separate structure function S m 2 (N, τ ) which are not necessarily independent because they are computed from the same data set. This protocol is similar to the jacknife and bootstrap methods employed in standard statistical analysis. The ζ m (N) we quote is the value obtained from an average over scaling exponents obtained for various discrete Eulerian point configurations.

Results
We first verify scalings for the second-order structure function at a single Eulerian point (N = 1) for different m-values in 3D. Figure 2(a) shows S m 2 (N = 1, τ ) for m = 1-4 as a function of the time τ normalized by the large eddy turnover time τ 0 ; all plots are vertically shifted for visual clarity. We find S m 2 (N = 1, τ ) ∼ τ 2/3 in accord with theory [17] and in agreement with published results [16,17,19] over nearly two decades in τ /τ 0 . As a confirmation, figure 2(b)  For comparison, we also plot the S m 2 (N = 1, τ ) to demonstrate the distinct difference in scalings, and the pre-factors for both functions are normalized to 1 to permit easy visual comparison. As is observed in figure 3, the structure functions scale as S m 2 (N → ∞, τ ) ∼ τ 2m/3 in the asymptotic limit. Finally, in figure 4 we plot the evolution of ζ m (N) versus N to track the convergence of the scaling exponent from point to the asymptotic field limit. We find the 3D convergence is exponential and obtain the best fit that follows the form where N 0 is the only fit parameter and represents the exponential convergence rate constant. We find N 0 varies within a narrow band between 24-29 points for all m-values reported here. Whereas we found no systematic dependence between m and the 24-29 point range, this range lies within the statistical variability of the various configurations over which the average exponent was calculated for the convergence of ζ m (N) versus number of points N. Before closing, we note in particular also captures the m = 1 case by default where the exponential vanishes and one trivially recovers ζ m (N) = 2/3 for all N values, and hence is consistent for all m.

Background
Turbulence in two spatial dimensions enjoys its own particularities in relation to its 3D counterpart. Fully developed 3D turbulence is characterized by a range of length scales, viz. the inertial range, where inertia strongly dominates over viscous effects, and indeed energy is a constant of motion within the inertial range in this inviscid limit (ν → 0). Kolmogorov's second seminal result of [2] tells us the inter-scale energy transfer within the inertial range of 3D turbulence proceeds or 'cascades' from large to small spatial scales. Below the Kolmogorov scale η, the energy is dissipated through viscous action. In contrast, 2D turbulence is characterized by the presence of two inviscid quadratic invariants, namely the mean kinetic energy and the mean enstrophy Ω or mean squared vorticity (Ω ≡ ω 2 , where ω ≡ ∇ × u is the vorticity) as was first elucidated by Kraichnan [30], Leith [31], and Batchelor [32] (KLB). The KLB theory explains that the conservation of mean kinetic energy and mean enstrophy leads to the presence of two simultaneous cascades in 2D turbulence. In this dual cascade scenario, the enstrophy cascades from the injection or forcing length scale r inj towards smaller scales; this is commonly referred to as the 'direct' or 'enstrophy' cascade regime. On the other hand, energy simultaneously cascades upscale from the injection scale to larger scales, and this regime is commonly called the 'inverse' cascade regime, to distinguish its direction from that of 3D turbulence. We refer the interested reader to several excellent reviews [33][34][35][36] and desist from further details except to state the primary scaling expectations for spectra in the direct and inverse cascade regimes. Also, we limit ourselves to velocity structure functions and stay away from vorticity structure functions since the current study concerns itself with fluctuations in velocity and its higher powers. The primary inertial range scaling expectations for velocity spectra in 2D turbulence are: where β ≡ d ω 2 dt , and around ω 2 denote a spatial average. The inverse cascade scaling of S 1 2 (r) ∼ r 2/3 is expected in the range of scales r inj < r < L, since energy is transferred from the injection scale r inj to larger scales until it spans the system size L, where it undergoes dissipation through large scale friction with the medium along the third spatial dimension [37]. Whereas equation (11a) is dimensionally consistent and has found both numerical and experimental validation [35,36], equation (11b) for the direct cascade regime is more involved. Specifically, S 1 2 (r) ∼ r 2 corresponds to E(k) ∼ k −3 spectrum which causes the enstrophy transfer rate to become logarithmically dependent upon the wavenumber k leading to an infrared divergence, i.e. as Re → ∞, the injection scale wavenumber k inj → 0 or r inj → ∞, which is unrealistic. Physically, this logarithmic dependence is attributed to the effect of larger structures contributing to the shear on smaller structures in the flow [36], thereby leading to non-local effects. Whereas a general theory rectifying this infrared divergence is missing, Kraichnan [38] proposed an approximation of the form E(k) ∼ β 2/3 k −3 [ln(k/k inj )] −1/3 such that the enstrophy transfer rate divergence remains suppressed. For this reason, a logarithmic correction factor ln(r inj /r) enters equation (11b), but experimentally difficult to measure due to the r 2 scaling to leading order.
Experimental measurements in soap films using both laser Doppler velocimetry (LDV), which represent a spatial point measurement in time that exploits TH, PIV representing an Eulerian field measurement, as well as PIV on electromagnetically forced thin fluid layers have all reported power-law scalings for S 2 (r) in the direct cascade regime. The scaling exponents reported in these measurements all fall in the range ζ 1 (N = 1) = 1.3-2.3 with the general consensus expectation, particularly from high-resolution numerical simulations, being about 1.6 [39][40][41][42]. Owing to several subtleties that arise from forcing as well as dissipation mechanisms in numerical schemes that are not relevant to our experiments, we refer the reader to [36] for details.
Moving beyond 2D turbulence spectra, Belmonte et al [43] experimentally verified the validity of TH in 2D for gravity driven decaying soap film turbulence. However to our knowledge, random sweeping effects have not been studied in 2D turbulence to date. Similarly, higher-order turbulence spectra have not been theoretically or experimentally explored in 2D turbulence to the best of our knowledge. In the absence of theor etical guidance and lacking in physical intuition, we resort to direct experimentation and proceed through comparative analysis with equivalent measurements in 3D.

Experimental methods
The 2D turbulence experiments were performed in a gravity driven soap film setup (figure 5) whose design closely followed that of prior studies [44][45][46][47]. The setup was comprised of two nylon fishing lines hung taut by a weight from a rigid metallic frame. Soap solution (2% mass concentration in MilliQ  deionized water) containing neutrally buoyant hollow glass tracer particles (pre-sieved to max 4 µm diameter, TSI Inc. Model 10089) for PIV imaging was injected onto the vertically hanging nylon wires through a nozzle. Although preliminary experiments were performed with a reservoir above the setup providing constant pressure head through a circulation system, we found by trial and error that a syringe pump injection method provided more reliable results in the cur rent study. The soap solution was injected using a 50 ml glass syringe with teflon cover pusher (to minimize stick-slip perturbations) mounted on a harvard apparatus (Elite II) syringe pump. The experiments were therefore conducted under constant flux rather than constant pressure, primarily to maintain film and flow stability and to avoid drift in flow conditions over the experimental duration. The nylon wires wetted by the draining soap solution were stretched apart at four points to generate a quasi-2D soap film of dimensions 2.10 m height, 4.13 cm width and based on prior studies, the thickness is of order 10 µm. A one dimensional grid (commercial hair comb) with 0.9 mm average tooth diameter and 0.82 ± 0.01 mm intertooth spacing was placed in the film perpendicular to the flow direction to generate turbulence in its wake. Turbulence intensity decays downstream with increasing distance from the grid, hence the measurement distance from the grid has a bearing on the results. All our measurements were performed in a window two channel widths below the grid to allow turbulence generated by the grid to fully develop. Prior works [36] report the turbulence becomes isotropic and homogeneous roughly one channel width below the grid. Air friction is expected to slow the flow to terminal velocity, but our measurements were made at a height where terminal velocity was not achieved.
The soap film was illuminated with two independent light sources for imaging purposes; one was a coherent laser sheet coplanar with the film (532 nm, 7.5 W coherent verdi diode pumped solid state laser) and a 20 W diffuse light source (Northstar 250W). A high speed camera (Vision Research Phantom v641) faced the soap film and image data was acquired at 5500 frames per second at a resolution of 800 × 800 pixels. The data was acquired for a square window of side length 2.2 cm in the film's central region leaving out roughly 0.9 cm on either side of the channel to avoid boundary effects. The raw image data was processed with DPIVSoft 2010 [48] to construct the velocity field with a grid spacing of 8 pixels. The spatial and temporal resolution in the 2D experiments are therefore 0.022 cm/pixel and 0.18 ms/frame respectively. The experiments were repeated under identical conditions eight times to collect sufficient temporal statistics. The relevant experimental parameters for the 2D experiments are listed in table 1. We note that the relevant length scale in the 2D experiments is the injection length scale r inj 1 mm, whereas the relevant time-scale in computation of structure functions is the flow correlation time τ 0 . The composite signal U m (t) was computed in the same manner for 2D and 3D data which we describe in section 2.1.
We close the experimental description with a few qualifiers. Firstly, although TH is known to apply in 2D [43], the turbulence intensity is high (please see table 1) in our experiments that we do not expect TH to hold for measurements discussed here. Secondly, it is well known that continuous forcing is necessary to study 2D turbulence in the inverse cascade regime [49] as is usually achieved with electromagnetically forced fluid layers [50], whereas decaying turbulence observed in gravity driven turbulence operates in the enstrophy or direct cascade regime [36]. All measurements we report concern velocity measurements in the direct cascade regime, a fact that has implications for interpretation of results to be discussed in the next section.
The computation of structure functions S m 2 (N, τ ), N 1 in 2D followed the same method detailed in section 2.2 with a few minor modifications. First, the mean velocity was subtracted from the 2D PIV velocity field, a step not required in 3D for lack of mean flow. However, considerable anisotropy still existed in the 2D PIV field from simple diagnostic tests performed for the streamwise and spanwise velocity components. The spherical harmonic decomposition to recover the isotropic sector from the PIV fields [29] necessary for our 3D data was also conducted on the 2D data. Unlike in 3D, for 2D measurements, the scaling exponent became shallower and approached a constant value at m = 5, measurements were taken up to m = 6. The rest of the procedure for computation of structure functions for both N = 1 and N > 1 cases is identical to the procedure followed for 3D data as detailed in section 2.2.

Results
The results for 2D direct cascade are very different from the 3D results presented above. In figure 6(a), we plot the 2D counterpart of figure 2, namely S m 2 (N = 1, τ ) for m = 1-4. First, the exponent we obtain for m = 1, S 1 2 (τ ) ∼ τ 1.54 is certainly below S 1 2 (τ ) ∼ τ 2 theory predicts but it is well within the range of ζ 1 (N = 1) = 1.3-2.3 reported by prior experiments [39][40][41][42] and numerical simulations [36]. We note that straightforward application of TH [39] is not possible in our measurements since I ≡ u rms /u ∼ 0.5, which in fact places our measurement within the RSH regime. Furthermore, with increasing m-value, we observe the scaling becomes progressively shallower in a systematic manner that is amenable to a quadratic fit of the form ζ m = 1.98 − 0.49m + 0.05m 2 (figure 6(c)) and reaches an asymptote around m = 5. The observed 2D behavior is clearly in sharp contrast with its 3D counterpart which yields ζ m (N = 1) = 2/3 for all m.
Next we present the asymptotic scalings where too the results diverge considerably from 3D behavior shown in figure 3. Figure 7 plots S m 2 (N, τ ) for m = 1-4 at N = 1 and N = 10 4 , the total number of PIV grid points and the respective log-derivatives. The scaling clearly steepens with increasing N. Looking at m = 1 scalings in figure 7(a), we observed ζ m (N) increases from ζ m (N = 1) = 1.54, close to prior observational reports [36,[39][40][41][42] and reaches ζ m (N tot. = 10 4 ) = 2.15 close to the theoretical prediction of τ 2 , ignoring logarithmic correction. To study the approach to asymptotic convergence, we plot ζ m (N) versus N in figure 8, this is the 2D counterpart for 3D result in figure 4. The best fit obtained is logarithmic of the form ζ m (N) = A + B · log(N) which improves by extending to smaller N values with increasing m, but starts deviating from log behavior and approaches a constant value beyond N = 1000. The logarithmic fit has two parameters A and B; ideally A = ζ m (N = 1) leaving only the convergence rate constant B as the lone fit parameter. However, we note the initially flat signature of ζ m (N) for the range 1 < N < 30 (varies with m value) in figure 8, and for this reason both A and B had to be treated as independent fit parameters for our plots. The measured values for logarithmic convergence are tabulated in table 2.

Discussion
The starting point for our discussion is the composite signal U m (t). Instantaneous, spatial fluctuations can capture the selfsimilar structure but say nothing about the nature of interscale energy transfer, which is a dynamical process and hence expressly necessitates a time-domain measurement. Since the sampled Eulerian points are spatially apart but concurrently sampled, the composite signal captures the inter-scale coupling for energy and/or enstrophy transfer through the composite signal's two-point correlation function: We see from equation (12) that the two-point correlator is decomposable into two terms. The first term on right hand side of equation (12) is the self-correlation term, it captures the correlation of u m (t) at ith Eulerian point with itself in time. The second term on the other hand encodes the spatial cross-correlation between temporal signals obtained from the ith and the j th locations. In other words, if the ith and j th Eulerian points are separated by a distance d, then temporal fluctuations over time scales corresponding to all eddies of size r < d decorrelate between u m i (t) and u m j (t), but all temporal fluctuations over time scales corresponding to eddies of size r d remain correlated. This picture is generally true in 3D and in the 2D inverse cascade regime, but becomes nuanced in the 2D direct cascade regime due to the non-local logarithmic correction factor discussed earlier. When d > l 0 all inter-scale correlations are lost and the cross-correlation term u m i (t)u m j (t + τ ) → 0, leaving behind only the self-correlation term. In the large N limit, the self-correlation term encodes the true temporal spectrum that is statistically independent from one scale to the other, and hence is expected to satisfy K41A expectations in 3D; as discussed earlier, no theoretical expectation is available for 2D.
For N = 1, the self-correlation term gives the temporal, Eulerian second-order structure function, there is no crosscorrelation term. But with increasing N, S m 2 (N, τ ) ∼ τ ζm(N) where ζ m (N) evolves from the Eulerian point towards a coarse-grained ( N 1) form of the field limit ( N → ∞). The evolution of ζ m (N) as a function of N is controlled by the cross-correlation term of equation (12) and it encodes the strength of the coupling between scales. In 3D turbulence, we expect inter-scale coupling to be local, in other words the energy transfer proceeds from length scale r in the inertial range to (r − δr) on average with some backscatter upscale towards (r + δr), and this inter-scale coupling strength falls exponentially. The reverse would hold true for the 2D inverse cascade regime where energy transfer proceeds towards large scales. The locality of inter-scale coupling is usually expressed in the wave number space, but holds true for real space as well. This locality of inter-scale coupling indeed forms the basis for scale filtered energy or enstrophy transfer analyses conducted in experiments [40,51] and large eddy simulations [52]. Due to exponential decay of inter-scale coupling, we expect the cross-correlation term in equation (12) to fall off exponentially in 3D and in the 2D inverse cascade regime. Consequently, we expect the evolution of ζ m (N) to be exponential in going from N = 1 → ∞.
The 2D direct cascade regime is not so straightforward on account of the fact that its spectrum, equation (11b) must be versus τ /τ 0 shows change in scaling with increase in N from N = 1 (solid symbols) to total number of points N tot. = 10 4 (empty symbols) for (a) m = 1 (black circles) with its log-derivative (b), (c) for m = 2 (blue squares) and (d) its logderivative, (e) m = 3 (red upright triangles) and (f) its log-derivative, and (g) m = 4 (green inverted triangles) with (h) its log-derivative. multiplied with the logarithmic correction factor ln(r inj /r). As discussed earlier, this logarithmic factor represents non-local inter-scale coupling due to large scale strains dominating over all smaller scales. As a consequence, the evolution of ζ m (N) cannot be expected to scale exponentially, but should be logarithmic. Be that as it may, the logarithmic dependence has never been directly observed in reported studies because the quadratic (r 2 ) scaling in equation (11b) always dominates over the logarithmic corrections, and that is to be expected.
Spatial spectra, be they in 2D or 3D can only demonstrate the self-similar structure of the eddies through their corresponding spatial fluctuations via instantaneous snapshots, they are not suitable to test for inter-scale couplings; for that a suitable combination of time and space spectra are required. One might argue the scale filtered enstrophy transfer [40] might have presented evidence of logarithmic inter-scale dependence but this expectation would be incorrect. The scale-filtering approach relies on applying a low pass filter to smooth all fluctuations corresponding to scales smaller than the filter scale. All fluctuations corresponding to fluctuations larger than the filter scale are passed through. If one focuses on a specific filter scale, it will certainly carry the logarithmic imprint of r inj but it cannot be decomposed out by filtering at a specific scale. Furthermore, if two filter scales were applied on the data to look at the flux past these two scales, it would subtract out the large-scale logarithmic factor. Ergo, the filterscale approach is not suited to test for logarithmic dependence in our opinion. On the other hand, the current approach we propose relies on the temporal spectral scaling exponent as a function of inter-scale separation, hence the approach of ζ m (N) could potentially exhibit logarithmic slow growth from N = 1 towards the asymptotic field limit N → ∞.
The 3D turbulence results presented above are in accord with existing theoretical understanding. For N = 1 case, and for all m-values, the results agree with prior theory and measurements [16,17,19], as are the asymptotic values of exponents which approach ζ m (N) = 2m/3 with increasing N. The exponential convergence we empirically determine from data  (table 2) is consistent with this picture given the non-local interaction the logarithmic correction factor introduces in the 2D direct cascade structure functions. Note that the conv ergence rate constant remains roughly constant in the range B = 0.4-0.5 for all m-values, whereas the N-value at which the logarithmic behavior commences exhibits monotonic decrease with m-value. This leads us to speculate that the convergence rate constant B may have physical origins, i.e. it may show dependence on Reynolds number or a similar variable. However the fit parameter A may point to statistical dependence.
We however caution the reader with a few qualifiers. Firstly, the logarithmic correction term has never been observed in experiments or numerics to our knowledge, and if indeed our interpretation of the second term in the composite correlation function of equation (12) is correct, further experimental and numerical evidence would be desirable to strengthen this interpretation. For instance, if similar measurements in inverse cascade regime in electromagnetically forced fluid layers do not exhibit exponential point to field convergence, our interpretation is rendered invalid. Also worthwhile exploring in this context are the vorticity structure functions in direct cascade regime, where Falkovich and Lebedev [53] predict a logarithmic behavior. Independent numerical evidence in both direct and inverse cascade regimes would go a long way in further probing this scenario since they have access to more quantities that are unavailable to experiments. However, it is difficult for simulations to simultaneously achieve high spatial and temporal resolution and when both are deemed necessary, numerical runs have been restricted to a duration spanning one turnover time. Finally, the logarithmic correction proposed by Kraichnan [38] comes about through imposition of the constant enstrophy flux requirement by fiat. A self-consistent field theoretic construction still eludes the community to this day,  for an extensive analysis and discussion we refer the reader to the excellent treatment by Bowman [54]. We therefore submit our interpretation is consistent with known behavior but not conclusive until verified in situations detailed above. The above qualifiers notwithstanding, some facts may still be gleaned. Firstly, the slow convergence reported in figure 8 cautions us as regards resolution requirements, both in experiments (PIV) and numerical simulations. Further work is needed on both low and large N limit behavior of the scaling exponent ζ m (N) in 2D. Given the soap film's channel boundaries in our experiment are roughly 0.9 cm from the imaging window edge, the large N deviation could possibly be a boundary effect and will be systematically explored in future. Secondly, as regards the relevance of these results to atmospheric flows and wind energy which formed our primary motivation in undertaking this study: atmospheric flows are subject to broadband forcing at scales larger than the small scales where wind energy is generated by turbines (maximum rotor diameter of order 100 m) and where wind speed measurements are conducted. Whereas our study is incapable of discriminating whether atmospheric flows are 2D, 3D, or mixed in character, at least at these small scales, we do expect the flows to fall within the direct cascade regime if atmospheric turbulence were to have 2D characteristics. The 3D (figure 4) and 2D (figure 8) convergence combined with the exponentially fast approach of scaling presented in [16] tells us at least at these small scales, it is safe to treat atmospheric boundary layer flows as three dimensional in character. Finally, an observation on the range ζ 1 for 2D. Figure 8 and table 2 show for 2D, ζ 1 (N) remains constant around 1.55 in the range N = 1-32 and approaches an asymptotic value close to 2 at large N. We note that values at the low end of the range ζ 1 = 1.3-2.3 [39][40][41][42] reported in literature come from LDV data representing N = 1 case whereas values at mid to high end of the range come from PIV data representing the field limit. We submit that this range of values for ζ 1 most probably arises from the discrepancy in Eulerian point versus field measurements, an issue that is absent by default in 3D.

Summary
We have presented experiments in 2D and 3D turbulence where we studied how higher-order spectra S m 2 (N, τ ) evolve from an Eulerian point (N = 1) towards a coarse-grained ( N 1) form of the asymptotic field limit N → ∞. Although the 3D measurements find consonance with prior theory and results and extend them further, the 2D results throw up several surprises. The logarithmic slow convergence in 2D as opposed to exponential convergence in 3D throws up an unexpected conundrum in particular. In spite of this, the results of this study do settle some of the questions arising from wind power fluctuations and atmospheric flows that formed the primary motivation for the study and will hopefully find use in the broader context of geophysical flows alike. In addition to serving a salutary warning for spatial resolution requirements in 2D experiments and simulations, the 2D results also shed light on why prior works report a wide range of exponents for velocity spectrum in the direct cascade regime. It is our earnest hope that these results will spur renewed interest in theoretical development, particularly in 2D turbulence.