Assessing Local Turbulence Strength from a Time Series

We study the possible link between “local turbulence strength” in a flow which is represented by a finite time series and a “chaotic invariant”, namely, the leading Lyaponuv exponent that characterizes this series. To validate a conjecture about this link, we analyze several time series of measurements taken by a plane flying at constant height in the upper troposphere. For each of these time series we estimate the leading Lyaponuv exponent which we then correlate with the structure constants for the temperature. In addition, we introduce a quantitative technique to educe the scale contents of the flow and a methodology to validate its spectrum.


Introduction
A deterministic definition of the atmospheric state as a function of time t in a domain D can be made in terms of the values of the various meteorological variables wind, temperature, pressure, moisture, etc. at each point of the domain.However, using this definition to characterize and make meaningful distinctions between different states is not practical due to lack of appropriate data.In view of this situation, one has to use averaged statistical quantities to characterize the atmospheric state at least partially e.g., one can use for this purpose the root-mean-square of any meteorological variable in D .Ideally a finite number of such quantities will be sufficient to characterize the state completely as in the case of an ideal gas in equilibrium .However it is obvious that this is not the case for the atmosphere.
In general, characterization of the state of a dynamical system in terms of a finite set of invariants is a difficult task.To complicate this problem further, the representation of the dynamical system is usually made in terms of a finite time series of measurements.For the atmosphere, the local state is usually represented by a finite time series of measurements at one point or almost simultaneous measurements over a spatial domain by aircraft .Under these circumstances, the challenge is to extract as much useful information from this representation 1, 2 .In particular one would like to use this time series to characterize and differentiate quantitatively between the different local atmospheric states.
One of the "traditional tools" for extracting information from this representation is the averaged spectrum of the time series and its slope at different wave numbers.This information is important for the determination of the atmospheric "structure constants" 3 , which impact the propagation of electromagnetic signals in the atmosphere 4 and many other applications.However, the estimation of this spectrum might be biased by discontinuities in the data and the parameters that are used for this estimation e.g., padding, windowing, etc. 5 .It is important therefore to cross-validate the determination of this spectrum by other independent methods.Our first task in this paper will be to use wavelet transform to carry out this cross validation.In addition, we introduce in this paper a new quantitative methodology to perform a "scale decomposition" of the flow from the time series data.
Our primary objective in this paper is to examine the possible characterization of the local turbulence strength 6 by a "Chaotic invariants" namely, the leading Lyaponuv exponent 1, 7, 8 of the time series that represents the flow.Thus, the basic conjecture that we want to validate in this paper, through multiple case study, is that local turbulence strength is functionally linked to the leading Lyaponuv exponent of the time series.Intuitively this means that as the turbulence strength in the flow increases, it becomes more chaotic and there is increased sensitivity to initial conditions.We provide some additional insights about this conjecture in Section 2.3.
To validate this conjecture we proceed indirectly.One well-known manifestation of turbulence strength in the atmosphere is through density fluctuations this leads to the well known "twinkling of the stars" phenomena .Hence turbulence strength is related to the value of C 2 N the refraction structure constant 3, 9, 10 .In the upper troposphere at heights of about 10 km , the main contributor to C 2 N is C 2 T -the structure constant for the temperature.It follows then that if our conjecture is correct then the leading Lyaponuv exponent for the temperature time series should be a function of C 2 T .We note that many other attempts were made to relate chaos theory and turbulence, see e.g., 11 .
In the analysis of atmospheric flow "length scale" is of great importance.In fact one speaks of the "integral scale" and the "scale regimes" that are present in the atmospheric flow.However, algorithms to compute these from a time series representation are available only for the integral scale 6 .Scale density representation which was introduced by Cohen 12 and others 13, 14 can give a quantitative representation of the different scales that are present in the flow and their intensity.In a way, this decomposition is similar to spectral analysis which represents the averaged energy density distribution as a function of frequency except that here we represent the intensity of the flow at a given scale.In particular, one can expect this intensity to "drop considerably" at the boundary between two scale regimes that are present in the flow.In this paper, we apply this transform to obtain a "scale decomposition" of the flow from its the time series data.
The plan of the paper is as follows: in Section 2 we present a short theoretical overview of the techniques that are used in this paper.In Section 3 we present the data that is being used to carry out the objectives of this research.In Section 4 we discuss the results that were obtained from this data using the techniques mentioned above.We end in Section 5 with some conclusions and observations.

Fractal Dimension and Wavelet Transform
In general, geophysical time series are assumed to be "self-affine" 5 .Accordingly one can apply the theorem of Mandelbrot and Van Ness which states that the Hausdorff dimension H d of the time series {y n } N 1 is related to the semivariogram by the relation where k is the "lag".Furthermore 1, 5 , the spectral density slope β of the time series is related to the Hausdorff and fractal correlation dimensions H d and D, resp.by the relation Together these relations enable us to evaluate β without recourse to the computation of the spectra.
Another method to compute β is based on the wavelet transform 5 where the wavelet function is chosen to be the Mexican hat function.Let W a, t be the wavelet transform of the time series y t where a is the "width" of the wavelet function.
It was found that for this wavelet the variance V a of W t, a satisfies the power law relation and furthermore which yields therefore another independent determination of β.
Finally we note that to evaluate β for a range of wave numbers where β is not constant a proper "detrending" namely, filtering has to be applied to the time series that is, one has to remove first the contributions from other parts of the flow.We discuss such a technique which is based on the principal component analysis in Section 2.4.

Scale Transform
The Fourier transform of a function f t is defined as This definition utilizes "weight functions" of the form e iωt which are eigenfunctions of the "frequency operator" Motivated by this observation Cohen 12 and De Sena and Rocchesso 13 introduced the symmetrized "scale operator" whose eigenfunctions are To interpret the parameter "c" in 2.10 , we note that s c, x has the same phase at positions x k , k 1, 2, . . .as x 0 if x k x 0 e 2πk/c .

2.11
The distance between two consecutive points in this sequence is not constant However, as "c" increases in value this distance shrinks; that is, s c, x oscillates more rapidly but with decreasing amplitude as x becomes larger .Thus "c" represents the "wave number" of the function s c, x which in this case measures the scale rather than frequency of this function.Observe that "meteorological larger scales" are represented by the smaller values of "c".Using these eigenfunctions, the scale transform of f x is defined as

2.13
It is easy to show using the transformation t e z that this integral is not singular at zero if f x is bounded.The scale density of f x is then defined as |D c | 2 .
To obtain some insight about the nature of this transform we, note that for an impulse Similarly for a wave

2.15
These are signals of "low" scale contents .However, for the eigenfunctions we have that is, these functions have a sharp scale contents.
In the context of the applications at hand, this transformation enables us to evaluate the contribution of different scales to a signal s x .

Structure Constants
The structure function of a geophysical variable, for example, the temperature T is defined as 3, 9, 10 where T are the turbulent fluctuations in the temperature and r is the vector from one point to another.Kolmogorov showed that for isotropic turbulence in the inertial range, this function depends only on d |r| and scales as C 2 T which appears as the proportionality constant in this equation is referred to as the "temperature structure constant".
The determination of the atmospheric structure constants 3, 9, 10, 15, 16 and in particular the temperature structure constant C 2  T is important in many applications, for example, the propagation of electromagnetic signals 3, 9 .Local peaks in the values of these constants, which are indicative of strong turbulence and reflect on the structure of the atmospheric flow, can have a negative effect on the operation of various optical instruments.
To estimate these structure constants in the upper troposphere or the stratosphere, it is a common practice to send high flying airplanes that collect data about the basic meteorological variables such as wind, temperature, and pressure along its flight path which may extend up to 200 kms.To estimate the averaged value of the structure constants along this path, one must decompose first the meteorological data into mean flow, waves, and turbulent residuals 17-19 .From the spectrum of the turbulent residuals, one can estimate an averaged value of the structure constants using Kolmogorov inertial range scaling and Taylor's frozen turbulence hypothesis 6 .For C 2 T in particular we have where F k is the temperature spectral density in the inertial range and k is the wave number.An averaged value for C 2 T over all wave numbers in the inertial range is obtained by averaging these values over k.Same relations hold for other geophysical variables.
In the upper troposphere where humidity is low, C 2 T is the main contributor to the refraction structure constant C 2 N 3, 9 which has an important impact on the operation of ground telescopes and is the cause for the twinkling of the stars.
To motivate our conjecture about the relationship between C 2 T and the leading Lyapunov exponents of the time series, we note that these exponents are defined by where f x 0 , t is a trajectory of the dynamical system with initial conditions x 0 and y is a small change in these conditions.Thus S r in 2.18 represents a spatial average of the function f r 1 r − f r 1 2 i.e., f x 0 , t T r, t .Observe that t was suppressed in 2.18 due to Taylor frozen turbulence field hypothesis .On the other hand, Lyaponuv exponents represent the asymptotic time behavior of the same function.The square is not material due the ln in the definition of the Lyaponuv exponents .Accordingly, our conjecture can be construed as an "ergodic proposition" which postulates a functional relationship between the spatial and time behavior of this function.

Data Decomposition
The statistical approach to turbulence splits the raw actual measurements of the flow variables u, T, p into a sum of

u u u u t , T T T T t , p p p p t , 2.22
where u, T, p represent the mean large scale flow; u , T , p represent waves are u t , T t , p t , "turbulent residuals" 17, 18 .
To effect such a decomposition in our data, we used the Karahunan-Loeve K-L decomposition algorithm or PCA which was used by many researchers.For a review see 20 .Here we will give only a brief overview of this algorithm within our context.Let be given a time series X of length N of some geophysical variable.We first determine a time delay Δ for which the points in the series are decorrelated.Using Δ we create n copies of the original series To create these one uses either periodicity or choose to consider shorter time series 21 .For the data under consideration, we used Δ 1024.For these n time series, one computes the Let λ 0 > λ 1 , . . ., > λ n−1 be the eigenvalues of R with their corresponding eigenvectors

2.25
The original time series T can be reconstructed then as where

2.27
The essence of the K-L decomposition is based on the recognition that if a large spectral gap exists after the first m 1 eigenvalues of R then one can reconstruct the mean flow or the large component of the data by using only the first m 1 eigenfunctions in 2.26 .A recent refinement of this procedure due to Penland et al. 20 is that the data corresponding to eigenvalues between m 1 1 and up to the point m 2 where they start to form a "continuum" represent waves.The location of m 2 can be ascertained further by applying the tests devised by Axford 22 and Dewan 16 .According to these tests turbulence data at the same location is a characterized by low coherence between u, v, w and a phase close to zero or π between w and T .A phase close to π/2 is characteristic of waves .These tests applied to our data show that to a large extent the residuals that were obtained from the K-L decomposition represent actual turbulence.For a detailed exposition of this decomposition with appropriate supporting plots see 18 .

Description of the Atmospheric Data
In this paper we will use two sets of data: the first was gathered by NASA 23 and the second by the US Airforce in collaboration with some Australian universities 10, 18, 19 .In Tables 1 and 2 we present the determinations of the spectral slopes which are based on 2.1 -2.6 for the raw data and the turbulent residuals.The estimate of the spectral slope based on fractal dimension and wavelet transform is between −2.5 and −3.On the other hand, the estimate based on the semivariogram is between −2. and −2.5.A possible explanation for this discrepancy is that these techniques give different weights to different parts of the spectrum namely, different wave numbers .Thus the semivariogram method gives more weight to higher wave numbers smaller scales where the slope should be −5/3 while the other two estimators give more weight to the larger scales where Bacmeister et al. found a slope of −3.Furthermore, the difference between the spectral slope of the different meteorological variable can be attributed to the fact that they are coupled differently to the stratospheric attractor in the sense of Lorenz 26 .We note also that the path of the plane taking the measurements traversed twice the polar vortex.

Scale Transform
We present in Figures 1 and 2 the scale densities for the temperature and w the vertical component of the wind for a flight that took place on September 9, 2002 over Australia at height of approximately 10 km above sea level.
These figures show clearly the integral scale s in the flow which is represented by small values of "c".They show also several other peaks which correspond to the smaller scales in the flow.The boundaries between the different scale regimes are represented by a well-defined dips in the scale density.We note also the similarities between the peaks of these two scale density representations for different values of "c".
The details that can be extracted from this analysis should be compared to the "traditional" statistical technique for estimating the integral scale only .This technique which is based on the zero crossing of the correlation function of the turbulent residuals is not always reliable as it might be sensitive to the methodology used to detrend the time series that represents the flow 6 .

C 2 T and the Leading Lyaponuv Exponent
To estimate C 2 T we used the US Airforce data.To derive the turbulent residuals, we used Karahunen-Loeve K-L decomposition as described in Section 2.4.For fifteen data sets that were contained in this data, averaged values for C 2 T were obtained using the methodology described in Section 2.3.
To compute the leading Lyaponuv exponents for these time series, we used Rosenstein algorithm and its implementation in the TISEAN package 7 .However, one should note that other algorithms are available for this purpose 7, 8 .To apply this algorithm, one has to determine first the "optimal" delay coordinates and embedding dimension.These were where E is the leading Lyaponuv exponents.Figure 3 presents a log-log plot of the results for this exponent versus C 2 T .We also show on this plot the least square line for this data.From this plot we see that a change in C 2 T over three orders of magnitudes correlates well with the leading Lyaponuv exponents.The fluctuations around the least squares line can be attributed to wave activity and possible measurements errors.This demonstrates that the Lyaponuv exponent can be used as a second local measure of turbulence strength in the data.In cases of discrepancy between C 2 T and the leading Lyaponuv exponent, one must trace out the reasons for this mismatch and correct them.

Summary and Conclusions
We introduced in this paper several data analysis tools that have the potential to validate and characterize the local atmospheric state and yield new insights about its structure.
In particular, we demonstrated that scale analysis can supplement other statistical and spectral tools which are used routinely in the analysis of meteorological data.We believe that we showed clearly the new depth that they bring into this analysis.
In addition, we showed that that there is a link between the local turbulence strength and the leading Lyaponuv exponent of the time series that represents the flow.Although our analysis does not constitute a proof of this conjecture in general, we still feel that our results will be useful in many applied contexts especially as a check for the validity of other global invariants that characterize the flow.Furthermore, the determination of the leading Lyaponuv exponent can help verify the value of C 2 T or other structure constants that have important practical applications.

Figure 1 :Figure 2 :
Figure 1: Scale transform for the temperature data for one of the flight segments on Sept 09, 2002 over Australia.

2 T 2 TFigure 3 : 2 T
Figure 3: Log-log plot of E-the leading Lyapunov exponent versus C 2 T .The first-order least squares fit for this data yield the relation E 10.85 * C 2 T 0.276 .

Table 1 :
Spectral slopes for raw data.