Mathematical Problems in Engineering

Volume 2010 (2010), Article ID 316841, 13 pages

http://dx.doi.org/10.1155/2010/316841

## Assessing Local Turbulence Strength from a Time Series

Department of Mathematical Sciences, Worcester Polytechnic Institute, 100 Institute Road, Worcester, MA 01609, USA

Received 20 January 2010; Accepted 24 March 2010

Academic Editor: Massimo Scalia

Copyright © 2010 Mayer Humi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

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.

#### 1. Introduction

A deterministic definition of the atmospheric state as a function of time in a domain 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 ). 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 (the refraction structure constant [3, 9, 10]). In the upper troposphere (at heights of about 10 km), the main contributor to is -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 . (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.

#### 2. Theoretical Background

##### 2.1. 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 of the time series is related to the semivariogram by the relation where is the “lag”. Furthermore [1, 5], the spectral density slope of the time series is related to the Hausdorff and fractal (correlation) dimensions ( and 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 ) be the wavelet transform of (the time series) where a is the “width” of the wavelet function.

It was found that for this wavelet the variance of 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.

##### 2.2. Scale Transform

The Fourier transform of a function is defined as This definition utilizes “weight functions” of the form 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 “” in (2.10), we note that has the same phase at positions as if
The distance between two consecutive points in this sequence is not constant
However, as “” increases in value this distance shrinks; that is, oscillates more rapidly but with decreasing amplitude (as becomes larger). Thus “” represents the “wave number” of the function 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 “”.

Using these eigenfunctions, the * scale transform* of is defined as
It is easy to show (using the transformation that this integral is not singular at zero if is bounded. The scale density of is then defined as .

To obtain some insight about the nature of this transform we, note that for an impulse Similarly for a wave (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 .

##### 2.3. Structure Constants

The structure function of a geophysical variable, for example, the temperature is defined as [3, 9, 10] where are the turbulent fluctuations in the temperature and is the vector from one point to another.

Kolmogorov showed that for isotropic turbulence in the inertial range, this function depends only on = and scales as

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 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 in particular we have where is the temperature spectral density in the inertial range and is the wave number. An averaged value for (over all wave numbers in the inertial range) is obtained by averaging these values over . Same relations hold for other geophysical variables.

In the upper troposphere where humidity is low, is the main contributor to the refraction structure constant [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 and the leading Lyapunov exponents of the time series, we note that these exponents are defined by where is a trajectory of the dynamical system with initial conditions and is a small change in these conditions. Thus in (2.18) represents a spatial average of the function (i.e., . Observe that 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 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.

##### 2.4. Data Decomposition

The statistical approach to turbulence splits the raw (=actual) measurements of the flow variables into a sum of where represent the mean (large scale) flow; represent waves are , “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 (of length ) of some geophysical variable. We first determine a time delay for which the points in the series are decorrelated. Using we create 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 . For these time series, one computes the auto covariance matrix Let be the eigenvalues of with their corresponding eigenvectors The original time series can be reconstructed then as where The essence of the K-L decomposition is based on the recognition that if a large spectral gap exists after the first eigenvalues of then one can reconstruct the mean flow (or the large component (of the data)) by using only the first eigenfunctions in (2.26). A recent refinement of this procedure due to Penland et al. [20] is that the data corresponding to eigenvalues between and up to the point where they start to form a “continuum” represent waves. The location of 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 and a phase close to zero or between and . (A phase close to 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].

#### 3. 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].

##### 3.1. NASA Data

This data was obtained during a mission to the arctic [23] by high flying ER2 plane over the northern polar vortex. Out this data we analyze in this paper two data sets. These sets represent the “state of the stratosphere” (traversed by the plane) on Feb 9 and 10, 1989. The choice of these two dates seems to be appropriate since according to mission records Feb 9 was a “normal day” while on Feb 10 the measurements were conducted at a time when “large mesoscale disturbances were observed”. The time series of measurements for each of these dates consists of over 100,000 data points for each of the geophysical variables and (which represent, respectively, the west-east, north-south, and vertical winds, temperature, and pressure). The measurements were taken at (almost) equal intervals of 0.2 sec (a spatial distance of 45 m).

##### 3.2. US Airforce Data

Previous to this data collection campaign, the high flying airplanes were equipped with only one meteorological probe. However, it was impossible to use this one-probe data to compute several important characteristics of the flow such as Brunt-Vaisala frequency and the Richardson number or to initiate simulations of the flow in the vicinity of the flight path (to determine the vorticity field and the scale structure of the flow). Furthermore, it was impossible to construct a dynamical model of the stratified flow and its related structure constants.

In order to overcome these shortcomings, a special purpose airplane was equipped with three probes (on the two wings and tail). This plane was used in the period of 1999–2002 to collect data over Australia and Japan about the geophysical flow at heights of about 8 km–12 km [10, 18, 19]. The data was collected at a frequency of 50 Hertz and spatial resolution of approximately 1 m.

Based on instrument specifications, the data noise should be at a relative error level of . This is confirmed by the eigenvalues obtained in the K-L decomposition where the last few eigenvalues (which reflect the noise level in the data) are of order of the leading eigenvalue.

#### 4. Results

##### 4.1. Determination of the Spectral Slope

From a meteorological (and fluid dynamics) point of view, a slope of in the “inertial range” of the spectrum [6] characterizes a three-dimensional turbulence. However, as was noted already by Lilly [24], this slope can change due to “geophysical factors”. Furthermore Kraichnan [25] showed that for 2-dimensional turbulence there is a part of the spectrum where this slope is Bacmeister et al. [15], who considered similar time series which represent the stratospheric flow, found using averaged spectral analysis that for a large range of wave numbers the spectral slope is close to Our analysis of the time series for the two dates under consideration tends to lend independent support for these findings. In fact our estimates for the spectral slopes are between and This implies that the stratospheric flow (which is strongly stratified) has some of the characteristics of 2D turbulence.

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 and On the other hand, the estimate based on the semivariogram is between and 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 while the other two estimators give more weight to the larger scales where Bacmeister et al. found a slope of 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.)

##### 4.2. Scale Transform

We present in Figures 1 and 2 the scale densities for the temperature and (the vertical component of the wind) for a flight that took place on September , 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 “”. 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 “”.

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].

##### 4.3. and the Leading Lyaponuv Exponent

To estimate 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 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 determined using the mutual information and false neighborhood algorithms [21] which are also implemented in the abovementioned package. This analysis led us to choose a four-dimensional embedding space with delay coordinates of data points. This led to the following (least squares) functional relationship: where is the leading Lyaponuv exponents. Figure 3 presents a log-log plot of the results for this exponent versus . We also show on this plot the least square line for this data. From this plot we see that a change in 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 and the leading Lyaponuv exponent, one must trace out the reasons for this mismatch and correct them.

#### 5. 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 (or other structure constants) that have important practical applications.

#### References

- H. D. I. Abarbanel, R. Brown, J. J. Sidorowich, and L. S. Tsimring, “The analysis of observed chaotic data in physical systems,”
*Reviews of Modern Physics*, vol. 65, no. 4, pp. 1331–1392, 1993. View at Publisher · View at Google Scholar · View at MathSciNet - J. M. Nese and J. A. Dutton, “Quantifying predictability variations in low order ocean-atmosphere model,”
*Journal of Climate*, vol. 6, pp. 185–204, 1993. View at Publisher · View at Google Scholar - R. L. Fanta, “Wave propagation in random media,”
*Progress in Optics*, vol. 22, pp. 342–399, 1985. View at Google Scholar - V. Tatarskii,
*Wave Propagation in Turbulent Medium*, Dover, New York, NY, USA, 1961. - B. D. Malamud and D. Turcotte, “Self-affine time series,”
*Advances in Geophysics*, vol. 40, pp. 1–90, 1999. View at Google Scholar - H. A. Panofsky and J. A. Dutton,
*Atmospheric Turbulence*, John Wiley & Sons, New York, NY, USA, 1989. - R. Hegger, H. Kantz, and T. Schreiber, “Practical implementation of nonlinear time series methods: the TISEAN package,”
*Chaos*, vol. 9, no. 2, pp. 413–435, 1999. View at Publisher · View at Google Scholar · View at PubMed · View at Zentralblatt MATH - M. T. Rosenstein, J. J. Collins, and C. J. De Luca, “A practical method for calculating largest Lyapunov exponents from small data sets,”
*Physica D*, vol. 65, no. 1-2, pp. 117–134, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - E. M. Dewan, “Optical turbulence forecasting,” Tech. Rep. AFGL-TR-80-0030, Air Force Geophysical Laboratory, 1980. View at Google Scholar
- M. Humi, “Estimation of the atmospheric structure constants from airplane data,”
*Journal of Atmospheric and Oceanic Technology*, vol. 21, no. 3, pp. 495–500, 2004. View at Publisher · View at Google Scholar - J. Peinke, M. Klein, A. Kittel, A. Okninsky, J. Parisi, and O. E. Roessler, “On chaos, fractals and turbulence,”
*Physica Scripta*, vol. T49, pp. 672–676, 1993. View at Google Scholar - L. Cohen, “The scale representation,”
*IEEE Transactions on Signal Processing*, vol. 41, pp. 3275–3293, 1993. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - A. De Sena and D. Rocchesso, “A fast mellin and scale transform,”
*EURASIP Journal on Advances in Signal Processing*, vol. 2007, Article ID 89170, 9 pages, 2007. View at Publisher · View at Google Scholar · View at Zentralblatt MATH - H.-F. Liu, Y.-Z. Yang, Z.-H. Dai, and Z.-H. Yu, “The largest Lyapunov exponent of chaotic dynamical system in scale space and its application,”
*Chaos*, vol. 13, no. 3, pp. 839–844, 2003. View at Publisher · View at Google Scholar · View at PubMed · View at Zentralblatt MATH · View at MathSciNet - J. T. Bacmeister, S. D. Eckermann, P. A. Newman et al., “Stratospheric horizontal wavenumber spectra of winds, potential temperature, and atmospheric tracers observed by high-altitude aircraft,”
*Journal of Geophysical Research D*, vol. 101, no. D5, pp. 9441–9470, 1996. View at Publisher · View at Google Scholar - E. M. Dewan, “On the nature of atmospheric waves and turbulence,”
*Radio Science*, vol. 20, no. 6, pp. 1301–1307, 1985. View at Publisher · View at Google Scholar - F. Einaudi and J. J. Finnigan, “Wave-turbulence dynamics in the stably stratified boundary layer,”
*Journal of the Atmospheric Sciences*, vol. 50, no. 13, pp. 1841–1864, 1993. View at Publisher · View at Google Scholar - M. Humi, “Approximate simulation of stratospheric flow from aircraft data,”
*International Journal for Numerical Methods in Fluids*, vol. 41, no. 2, pp. 209–223, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - T. J. Duck and J. A. Whiteway, “The spectrum of waves and turbulence at the tropopause,”
*Geophysical Research Letters*, vol. 32, no. 7, Article ID L07801, 2005. View at Publisher · View at Google Scholar - C. Penland, M. Ghil, and K. M. Weickmann, “Adaptive filtering and maximum entropy spectra with application to changes in atmospheric angular momentum,”
*Journal of Geophysical Research*, vol. 96, no. D12, pp. 22659–22671, 1991. View at Google Scholar - M. B. Kennel, R. Brown, and H. D. I. Abarbanel, “Determining embedding dimension for phase-space reconstruction using a geometrical construction,”
*Physical Review A*, vol. 45, no. 6, pp. 3403–3411, 1992. View at Publisher · View at Google Scholar - D. N. Axford, “Spectral analysis of an aircraft observation of gravity waves,”
*Quarterly Journal of the Royal Meteorological Society*, vol. 97, pp. 313–321, 1971. View at Publisher · View at Google Scholar - S. Hipskind and S. Gains, “Airborne Arctic Stratospheric Expedition, NASA/VARP003 Appendix A: ER2 Meterorological Measurement Supplement 5Hz Data,” 5530-CDROM-283000, 1991.
- D. K. Lilly, “Numerical simulation of two-dimensional turbulence,”
*Physics of Fluids*, vol. 12, no. 12, pp. 240–249, 1969. View at Google Scholar · View at Zentralblatt MATH - R. H. Kraichnan, “Inertial ranges in two-dimensional turbulence,”
*Physics of Fluids*, vol. 10, no. 7, pp. 1417–1423, 1967. View at Publisher · View at Google Scholar - E. N. Lorenz, “Dimension of weather and climate attractors,”
*Nature*, vol. 353, no. 6341, pp. 241–244, 1991. View at Publisher · View at Google Scholar