Assessing Local Turbulence Strength from a Time Series

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𝛾_𝑘=1𝑁_𝑛(𝑦(𝑛+𝑘)−𝑦(𝑛))² (2.1) by the relation𝛾_𝑘=𝑘^𝐻𝑑, (2.2) 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−𝛽=2𝐻_𝑑+1=5−𝐷. (2.3)

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) 𝑦(𝑡)1𝑊(𝑡,𝑎)=2𝜋𝑎^∞_−∞𝑔(𝑎,𝑡−𝑡^′)𝑦(𝑡^′)𝑑𝑡^′,1𝑔(𝑎,𝑡)=()2𝜋^1/2(1−𝑥²)𝑒^−𝑥2/2𝑡,𝑥=𝑎, (2.4) where a is the “width” of the wavelet function.

It was found that for this wavelet the variance 𝑉_𝑎 of 𝑊(𝑡,𝑎) satisfies the power law relation𝑉_𝑎=𝑎^𝐻𝑤 (2.5) and furthermore𝛽=−𝐻_𝑤 (2.6) 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 as1𝑓(𝜔)=√2𝜋^∞−∞𝑓(𝑡)𝑒^{−𝑖𝜔𝑡}𝑑𝑡. (2.7) This definition utilizes “weight functions” of the form 𝑒^𝑖𝜔𝑡 which are eigenfunctions of the “frequency operator”1ℱ=𝑖𝑑.𝑑𝑡 (2.8) Motivated by this observation Cohen [12] and De Sena and Rocchesso [13] introduced the symmetrized “scale operator” 𝒞_𝑥=1𝑑2𝑖(𝑥+𝑑𝑑𝑥𝑑𝑥𝑥) (2.9) whose eigenfunctions are1𝑠(𝑐,𝑥)=√𝑒2𝜋√𝑥^{𝑖𝑐ln𝑥}. (2.10)

To interpret the parameter “𝑐” in (2.10), we note that 𝑠(𝑐,𝑥) has the same phase at positions 𝑥_𝑘,𝑘=1,2,… as 𝑥₀ if 𝑥_𝑘=𝑥₀𝑒^{2𝜋𝑘/𝑐}. (2.11) The distance between two consecutive points in this sequence is not constant 𝑥_𝑘+1−𝑥_𝑘=𝑥₀(𝑒^2𝜋/𝑐−1)𝑒^{2𝜋𝑘/𝑐}. (2.12) 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 1𝐷(𝑐)=√2𝜋^∞0𝑒𝑓(𝑥)^{−𝑖𝑐ln𝑥}√𝑥𝑑𝑥. (2.13) 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𝑓₁(𝑥)=𝛿(𝑥−𝑥₀1),𝐷(𝑐)=√𝑒2𝜋√𝑥^{𝑖𝑐ln𝑥0}₀. (2.14) Similarly for a wave𝑓₂(𝑥)=𝑒^{−𝑖𝜔0𝑥},𝑓(𝜔)=𝛿(𝜔−𝜔₀1),𝐷(𝑐)=√𝑒2𝜋√𝜔^{−𝑖𝑐ln𝜔0}₀. (2.15) (These are signals of “low” scale contents). However, for the eigenfunctions𝑠(𝑥)=exp(𝑖𝑐₀ln𝑥)√,2𝜋𝑥 (2.16) we have𝐷(𝑐)=𝛿(𝑐−𝑐₀), (2.17) 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]𝑆(𝐫)=⟨[𝑇^′(𝐫₁+𝐫)−𝑇^′(𝐫₁)]²⟩, (2.18) 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𝑆(𝑑)=𝐶²_𝑇𝑑^2/3. (2.19)

𝐶²_𝑇 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𝐶²_𝑇(𝑘)=4𝐹(𝑘)𝑘^5/3, (2.20) 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𝛼_𝑖=lim_𝑡→∞1𝑡ln‖𝐟(𝐱₀+𝐲,𝑡)−𝐟(𝐱₀,𝑡)‖,‖𝐲‖ (2.21) 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 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.

2.4. Data Decomposition

The statistical approach to turbulence splits the raw (=actual) measurements of the flow variables 𝐮,𝑇,𝑝 into a sum of𝐮=𝐮+𝐮^′+𝐮_𝑡,𝑇=𝑇+𝑇^′+𝑇_𝑡,𝑝=𝑝+𝑝^′+𝑝_𝑡, (2.22) 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 𝑋(𝑘),𝑋(𝑘+Δ),…,𝑋(𝑘+(𝑛−1)Δ). (2.23) 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 𝑛 time series, one computes the auto covariance matrix 𝑅=(𝑅_𝑖𝑗)𝑅_𝑖𝑗=^𝑁_𝑘=1𝑋(𝑘+𝑖Δ)𝑋(𝑘+𝑗Δ). (2.24) Let 𝜆₀>𝜆₁,…,>𝜆_𝑛−1 be the eigenvalues of 𝑅 with their corresponding eigenvectors 𝜙^𝑖=(𝜙^𝑖₀,…,𝜙^𝑖_𝑛−1),𝑖=0,…,𝑛−1. (2.25) The original time series 𝑇 can be reconstructed then as𝑋(𝑗)=^𝑛−1_𝑘=0𝑎_𝑘(𝑗)𝜙^𝑘₀, (2.26) where𝑎_𝑘1(𝑗)=𝑛^𝑛−1_𝑖=0𝑋(𝑗+𝑖Δ)𝜙^𝑘_𝑖. (2.27) 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 𝑚₁+1 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 𝜋/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].

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 10⁻³. 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 10⁻³ 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 −5/3 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 −3. 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 −3. 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 −2 and −3. 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 −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.)

Table 1 

Spectral slopes for raw data.

Table 2 

Spectral slopes for turbulent residuals.

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 9, 2002 over Australia at height of approximately 10 km above sea level.

Figure 1 

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

Figure 2 

Scale transform for $w$ (vertical wind) data for the same flight segment as in Figure 1.

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 1024 data points. This led to the following (least squares) functional relationship: 𝐸=10.85(𝐶²_𝑇)^0.276, (4.1) 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.

Figure 3 

Log-log plot of $E$—the leading Lyapunov exponent versus $C_T^2$. The first-order least squares fit for this data yield the relation $E=10.85*{(C_T^2)}^{0.276}$.

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.