## Abstract

The climatology of the spatial structure functions of velocity and temperature for various altitudes (pressure levels) and latitude bands is constructed from the global rawinsonde network and from Aircraft Communications, Addressing, and Reporting System/Aircraft Meteorological Data Relay (ACARS/AMDAR) data for the tropics and Northern Hemisphere. The ACARS/AMDAR data provide very dense coverage of winds and temperature over common commercial aircraft flight tracks and allow computation of structure functions to scales approaching 1 km, while the inclusion of rawinsonde data provides information on larger scales approaching 10 000 km. When taken together these data extend coverage of the spatial statistics of the atmosphere from previous studies to include larger geographic regions, lower altitudes, and a wider range of spatial scales. Simple empirical fits are used to approximate the structure function behavior as a function of altitude and latitude in the Northern Hemisphere. Results produced for spatial scales less than ∼2000 km are consistent with previous studies using other data sources. Estimates of the vertical and global horizontal structure of turbulence in terms of eddy dissipation rate *ε* and thermal structure constant *C _{T}*

^{2}are derived from the structure function levels at the smaller scales.

## 1. Introduction

Climate change studies have traditionally concentrated on the mean statistical properties of atmospheric variables. Early work was based primarily on rawinsonde data (Lorenz 1967; Palmen and Newton 1969) while more recent work uses global numerical weather prediction (NWP) model reanalyses (Hoskins et al. 1989; Kalnay et al. 1996; Uppala et al. 2005) to evaluate long-term global mean statistics (Fil and Dubus 2005; Peters et al. 2008). But other statistical properties of global atmospheric variables, such as long-term energy spectra and spatial correlations, are also of interest. For instance, climatologies of the observed spatial statistics are essential for rigorous evaluation of the effective spatial resolution of NWP models (Frehlich and Sharman 2004, 2008), for correct interpretation of forecast error statistics (Frehlich 2008), and for optimal data assimilation (Frehlich 2006).

In the seminal paper by Nastrom and Gage (1985) it was shown that the average spatial spectrum of atmospheric kinetic energy and temperature obeys a *k*^{−5/3} behavior from scales of a few kilometers to hundreds of kilometers, transitioning to a roughly *k*^{−3} behavior at longer wavelengths (*k* is the one-dimensional horizontal wavenumber). These results were based on measurements from specially instrumented commercial aircraft from the Global Atmospheric Sampling Program (GASP). More recent data from the Measurements of Ozone and Water Vapor by In-Service Airbus Aircraft (MOZAIC) campaign confirmed this behavior (Lindborg 1999; Cho and Lindborg 2001) and it has also been duplicated with high-resolution general circulation models (GCMs) and mesoscale NWP models (Koshyk and Hamilton 2001; Skamarock 2004; Frehlich and Sharman 2004; Takahashi et al. 2006; Hamilton et al. 2008).

However, the true climatology of the spatial statistics remains uncertain, since the inference from these statistics is somewhat incomplete. The GASP and MOZAIC sampling is limited to commercial aircraft flight tracks and cruising altitudes (9–12 km), and since high turbulent regions are intentionally avoided, the statistics inferred from these observations may be biased low compared to the true climatology. This problem can be avoided using NWP reanalysis data, however, the models suffer from changes in available data and data assimilation techniques (Bengtsson et al. 2004), and from model filtering effects that tend to underrepresent the smallest scales (Skamarock 2004; Frehlich and Sharman 2008).

Here we produce a more complete climatology based on data derived from a combination of global sounding data (SND) for the larger scales and routine automated winds and temperature measurements from a variety of commercial aircraft [Aircraft Communications, Addressing, and Reporting System (ACARS) and the Aircraft Meteorological Data Relay (AMDAR) system] to improve the climatological estimates at the smaller scales. This avoids the problems associated with the use of model data and extends the sampling regions provided by the GASP and MOZAIC data to higher and lower latitudes and altitudes.

Here, the derived climatologies of spatial variability are based on the average shape and level of structure functions computed from pairs of observations of temperature and velocity from rawinsonde and ACARS/AMDAR measurements for various altitudes, latitude bands, seasons, and geographic locations over the Northern Hemisphere. Structure functions are better suited for the analysis of spatial statistics given the unevenly spaced observations available and can easily be converted to spatial spectra (Frehlich and Sharman 2008). Structure functions were used by Lindborg (1999) and Cho and Lindborg (2001), and have been used for other meteorological applications (Buell 1960; Barnes and Lilly 1975; Maddox and Vonder Haar 1979; Gomis and Alonso 1988). Structure functions also permit local estimates of small-scale turbulence (Frehlich and Sharman 2004) based on recent theoretical interpretations of stably stratified turbulence (Lindborg 1999, 2006; Riley and Lindborg 2008).

The remainder of the paper is organized as follows: section 2 describes the data used. Section 3 describes the data analysis techniques and develops empirical fits to the structure functions. Section 4 presents results for different altitudes, regions, and seasons. Section 5 concludes the paper with a discussion and summary and implications for current theories of mesoscale turbulence. However, we do not attempt to explain the observed scalings produced here. This is an active area of research (Cho et al. 1999; Lindborg 1999, 2006; Lindborg and Cho 2001; Tung and Orlando 2003; Tulloch and Smith 2006, 2009; Gkioulekas and Tung 2006; Riley and Lindborg 2008) that is beyond the scope of this paper.

## 2. Data sources

The rawinsonde network is concentrated over continental areas of the Northern Hemisphere (see Fig. 1.4.2 of Kalnay 2003) where the average spacing between stations is ∼300 km. Here we use an archive (http://dss.ucar.edu/datasets/ds353.4/) of all rawinsonde measurements taken over 34 yr (1973–2006) for the standard times of 0000 and 1200 UTC. The measurement errors are ∼0.5 K for temperature (Luers 1997) and ∼0.5 m s^{−1} in each horizontal velocity component (Jaatinen and Elms 2000).

The ACARS and AMDAR observations of wind and temperature are used for the time period from 1 January 2001 to 31 December 2008. Here we use the term “AMDAR” for both of these datasets collectively. The sampling rate varies from as little as 10 s to as much as 30 min and the accuracy of each horizontal wind component is *σ _{υ}* ≈ 1.25–1.5 m s

^{−1}and of temperature is

*σ*≈ 0.5 K (Benjamin et al. 1999; Drüe et al. 2008). The discretization of the AMDAR data is 0.5 m s

_{T}^{−1}for wind speed, 1° for wind direction, and 0.1–0.2 K for temperature (Benjamin et al. 1999) and therefore the contribution of random discretization error to total instrument error is negligible. All the AMDAR analyses presented here used the Meteorological Assimilation Data Ingest System (MADIS)-formatted archive (see http://madis.noaa.gov/ for more details). Although more extensive than the GASP or MOZAIC datasets, the AMDAR data are still limited to typical commercial aircraft flight paths; however, the coverage is nearly global, and contains measurements taken during takeoffs and landings, and thus should provide sufficient coverage to estimate the climatology as a function of location, altitude, and time of year. Examples of the AMDAR data locations for one day at two different aircraft flight levels (FL; isobaric surfaces) are shown in Fig. 1. The 250-hPa region (224–276 hPa, which corresponds to commercial aircraft flight levels of approximately 32 000–36 000 ft or 9.75–11.0 km MSL) has the highest concentration of data (Fig. 1a). But even the smaller volume of data provided by AMDAR-equipped aircraft during climbs and descents (e.g., the 500-hPa region corresponding to an FL of approximately 18 000 ft or 5.5 km MSL shown in Fig. 1b) provides ample observations for robust determination of spatial statistics when the entire 7-yr observation period is used (from 109 000 measurements per day in 2002 to 212 000 measurements per day in 2007). However, as can be seen in Fig. 1, the AMDAR data are densest over the United States and Europe and therefore, like the rawinsonde network, do not provide a uniform sampling over the globe. Nevertheless, there is still enough coverage over the North Atlantic and North Pacific to construct fairly robust climatologies over most of the Northern Hemisphere.

## 3. Structure function estimates

For any pair of observation points of variable *q* separated by a horizontal vector **s**, the structure function *D _{q}*(

**s**) is given by {e.g., Monin and Yaglom [1975, p. 94, Eq. (13.40)]}

where 〈 〉 denotes an ensemble average. Figure 2 shows the geometry of the structure function analysis for two-point observations of velocity with east–west velocity component *u* and north–south velocity component *υ*. Two point observations of temperature share the same geometry. The spatial coordinates of each observation *i* are given by its spherical coordinates (*ϕ _{i}*,

*λ*,

_{i}*R*), where

_{i}*ϕ*,

_{i}*λ*, and

_{i}*R*denote latitude, longitude, and distance from the center of the earth, respectively. Latitude and longitude are expressed in radians and

_{i}*R*=

_{i}*R*+

_{E}*h*, where

_{i}*R*is the radius of the earth and

_{E}*h*is the altitude (MSL) of the observation.

_{i}For rawinsonde data, the latitude and longitude are adjusted for balloon drift using a constant rise rate of 5.0 m s^{−1} (Luers 1997; Luers and Eskridge 1998) and the recorded vector winds. Only rawinsonde data with complete wind speed profiles are used in the analysis to ensure reliable drift calculations. For ease of computation, spatial statistics are calculated on isobaric surfaces defined by the mandatory pressure levels provided in the sounding archive (700, 500, 400, 300, 250, and 200 hPa). The statistics derived for constant pressure surfaces are equivalent to those derived for constant altitude levels (Skamarock 2004).

For each pressure level, the difference between the spatial coordinates of any two observations on the sphere with mean distance from the center of the earth *R* = (*R*_{1} + *R*_{2})/2 is given by the arc distance vector **s** = (*s _{x}*,

*s*) with components defined by (Fig. 2)

_{y}where Δ*λ* = *λ*_{2} − *λ*_{1}, Δ*ϕ* = *ϕ*_{2} − *ϕ*_{1}, and *ϕ _{a}* = (

*ϕ*

_{1}+

*ϕ*

_{2})/2. The arc distance

*s*is along constant latitude with positive to the east and

_{x}*s*is arc distance along constant longitude with positive to the north. The heading

_{y}*ψ*of observation 2 with respect to observation 1 is

Structure functions of the velocity field are calculated for the longitudinal component *υ _{L}* (along the separation vector) and the normal component

*υ*(normal or transverse to the separation vector), that is,

_{N}where *u _{i}* and

*υ*are the east and north velocity components and

_{i}*p*denotes the pressure level. This convention was used by Buell (1960) and other early studies, but other authors have chosen different separation variables and different conventions for the transverse velocity orientation (Lindborg 1999; Cho and Lindborg 2001). However, the structure functions’ estimates were found to depend only weakly on the convention used. For example, the difference between this formulation and that of Lindborg (1999) at a separation of 1000 km is less than 3% and at a separation of 400 km is less than 1%.

To ensure a sufficient number of samples for stable results, the velocity and temperature structure functions are computed as averages over longitude bands between *λ*_{min} and *λ*_{max} (usually 0°–360°) and within latitude bands *ϕ*_{min}–*ϕ*_{max} (usually with a bandwidth of 10°). That is, for the velocity structure functions, Eq. (1) is implemented as

where *q* is either *L* or *N*, and *D _{LL}* and

*D*will denote velocity structure functions throughout,

_{NN}*s*= (

*s*

_{x}^{2}+

*s*

_{y}^{2})

^{1/2}is the arc distance between the two observations, and

*M*denotes the total number of observation pairs in the latitude band. The separation bin around

*s*is of a decade in log space. However, the mean temperature

*T*(

*p*,

*ϕ*,

*λ*) can have a pronounced variation with latitude away from the tropics and must be removed. Therefore, the estimate for the temperature structure function

*D*subtracts out the meridionally averaged mean over the latitude band:

_{T}where

For *D _{LL}* and

*D*, the 0000 and 1200 UTC sounding data from 1973 to 2006 are used. Since the temperature sensor from rawinsondes changed after 1988 (Luers and Eskridge 1998), only the time period 1989–2006 is used for the temperature structure functions. The sounding data sample the larger scales and therefore the contribution from the observation error is small compared with the large temperature differences in the atmosphere.

_{NN}To improve the statistical accuracy of the structure functions from AMDAR data, especially for lower altitudes, each altitude region was divided into 13 pressure intervals of depth 4 hPa for a total extent of 52 hPa centered on each of the mandatory pressure levels of the sounding data (same processing domain as in Fig. 1). The structure functions were calculated for each 4-hPa pressure interval and then averaged to produce a final estimate. Although we assume there are no effects from time differences between soundings taken at different locations, the time difference between the AMDAR observation pairs is restricted to be less than 1 h. This effectively limits the largest separations that can be computed using AMDAR data to ∼1000 km. Only data from the *same* aircraft are processed to reduce the contribution of the random bias from different aircraft (Ballish and Kumar 2008; Drüe et al. 2008). This is consistent with the approach used by Lindborg (1999) and Cho and Lindborg (2001), for the GASP and MOZAIC data.

The average structure functions are produced using only sounding data with high-quality flags (“A” or “space”). Similarly, only AMDAR data that passed all quality control (QC) stages are used (Moninger et al. 2003). However, examination of the quality controlled AMDAR data showed some remaining inconsistencies between the recorded times and locations for some observation pairs. Therefore, structure functions are produced from AMDAR data pairs only when the effective ground speed, *s*/Δ*t*, is between 50 and 400 m s^{−1}. This criterion is based on consideration of typical aircraft airspeeds of 100–250 m s^{−1} and expected maximum head (tail) winds of 50 (150) m s^{−1}.

Examples of the probability density function (PDF) from the rawinsonde and AMDAR data are provided in Figs. 3 and 4. Figure 3 shows the PDF of wind speed and temperature derived from the sounding data. Outliers in the temperature data are clearly evident around the dominant atmospheric feature. There is a small feature above 100 m s^{−1} in the wind speed PDFs. These outliers are removed using a threshold. We also assume that the impact of any remaining bad data (the approximately constant floor of the PDF around the peak) is negligible. This is a good approximation if the total fraction of bad data accepted by the threshold level is small.

For small separations, the structure functions have small values and are sensitive to measurement errors, especially for the AMDAR data. This is demonstrated in Fig. 4 where the tails of the PDF of the velocity differences squared of the AMDAR data appear to deviate from the local scaling in the tails of the distribution. Because of the poor statistical accuracy produced by the small number of events in the extreme tails of the PDF, it is difficult to distinguish legitimate rare events with large horizontal velocity or temperature gradients from random outliers. Thus, in addition to the location and time consistency checks mentioned above, a maximum threshold of 1.5 m s^{−1} km^{−1} was applied to the horizontal wind shear estimates Δ*u*/*s* and Δ*υ*/*s* to remove apparently bad wind data. Most of these cases were due to a known error in the AMDAR data wind direction noted by Pauley (2002), and this algorithm successfully removes these errors. Similarly, data with horizontal temperature gradients Δ*T*/*s* > 0.4 K km^{−1} were also removed.

Examples of global velocity and temperature structure functions for the AMDAR and sounding data are shown, respectively, in Figs. 5 and 6 for the 250-hPa pressure level at 40°–50°N. Both the velocity and temperature structure functions exhibit the expected *s*^{2/3} behavior at small lags (equivalent to the *k*^{−5/3} scaling of spectra at high wavenumber *k*). The AMDAR temperature data typically intersect the instrument noise floor for spacings less than ≈8 km (not shown), so a minimum separation of 10 km for temperature is used (Fig. 6). However, the velocity data appear to be less noisy at the small scales, and reasonable structure functions are obtained for lags as small as 3 km (Fig. 5). The estimation error is largest for the small lags since they have the smallest number of data pairs; this is clearly shown by the scatter in the sounding structure functions at ∼100-km spacing and in the AMDAR velocity structure functions at lags smaller than ∼3 km. Nevertheless, the agreement of the AMDAR and sounding structure functions in the overlap regions of ∼100–400 km is surprisingly good considering the different spatial and temporal averaging used (see Fig. 1; Kalnay 2003, her Fig. 1.4.2).

For comparison, the Lindborg (1999) average model is also shown in Fig. 5. Note that there is no corresponding model for comparison to the temperature structure functions shown in Fig. 6. The shape of the velocity structure functions is quite consistent with the Lindborg model at all scales, although the AMDAR measurements give levels slightly lower than Lindborg’s at the smaller scales. The level of the transverse structure function *D _{NN}* is higher than the level of the longitudinal structure function

*D*at the larger lags (50 km <

_{LL}*s*< 2000 km), mainly because of the sensitivity of

*D*at these larger lags to meridional wind variations induced by planetary waves. Lindborg (1999) has also shown that

_{NN}*D*at lags between 200 and 1000 km is in good agreement with the predictions of 2D isotropic turbulence based on

_{NN}*D*. At smaller lags, the

_{LL}*D*and

_{NN}*D*merge and consequently, as discussed by Lindborg (1999), neither the 2D isotropic scaling (Ogura 1952; Lindborg 1999) nor the 3D isotropic scaling (Monin and Yaglom 1975, p. 353) are valid.

_{LL}The ground speed and velocity shear thresholds applied to the AMDAR data were tested by varying the thresholds until little change was produced in the resulting structure functions and the *s*^{2/3} scaling was consistently observed at small lags. We therefore believe the structure functions presented here are dominated by atmospheric processes and not observation errors, which would produce a constant value at small lags. From the values of the velocity structure function at spacings less than ∼2 km in low-turbulence regions, the AMDAR random measurement error is estimated to be less than 0.4 m s^{−1} (see Fig. 5 at small lags). The total measurement error based on differences between two different AMDAR aircraft is difficult to estimate because of the atmospheric contribution. However, we believe the total AMDAR observation error of 1.25–1.5 m s^{−1} per velocity component estimated by Benjamin et al. (1999) and Drüe et al. (2008) may be biased high because of the contribution of atmospheric turbulence.

As shown in Figs. 5 and 6, the structure functions have a simple shape and therefore can be fit to simple empirical models as demonstrated by Lindborg (1999). We use the following empirical model that includes a better representation of the larger scales:

where *k* is a wave modal number, *J* is the number of modes, and mod is *LL*, *NN*, or *T*. The coefficients *a _{k}* describe the smaller-scale turbulence while the random amplitudes of the larger-scale planetary waves are modeled by the cosine terms with amplitude

*b*and wavelength

_{k}*c*. This empirical model satisfies the ubiquitous

_{k}*s*

^{2/3}scaling for small lags and includes the larger-scale planetary waves and baroclinic waves. The deviation from the

*s*

^{2/3}scaling to an approximate

*s*

^{a3}scaling occurs at

*a*

_{2}and for most cases,

*a*

_{2}<

*a*

_{4}, where

*a*

_{4}defines the largest turbulent scale before planetary waves become important. The tropics and some of the lower altitude regions appear to have a single-scale model (i.e., no planetary-wave component).

The best estimates of the average structure function are produced by combining the information from both the AMDAR and sounding data. We assume that the sounding data provide the best estimate of the larger scales while the AMDAR data provide more accurate small-scale statistics. The best-fit model is therefore produced by minimizing the error

where *a _{i}* denotes the parameters of the empirical model in Eqs. (10)–(11),

*J*is the number of AMDAR estimates, and

_{A}*J*is the number of sounding estimates. The estimation error of the structure function estimates is proportional to the mean since it is dominated by the turbulent processes (Lenschow et al. 1994). Therefore, the log of the structure function estimates has approximately constant error. Structure function estimates with high estimation errors resulting from small sample sizes are removed by using a threshold based on the relative error, typically 0.04 for the AMDAR data and 0.02 for the sounding data. The Powell algorithm (Press et al. 1986) produces the minimization.

_{S}The best-fit models and the structure-function estimates that passed the threshold test are shown in Figs. 7 and 8 for the same region as shown in Figs. 5 and 6, respectively, as well as a lower altitude of 700 hPa. Although these results are for a limited altitude and latitude region, the velocity structure functions agree well with the Lindborg (1999) model. The AMDAR data resolve the *s*^{2/3} scaling at small scales and the results indicate that the *s*^{2/3} scaling is valid over a large altitude and spatial region, especially for *D _{LL}*. The higher level of

*D*compared with

_{NN}*D*at 1000-km scales reflects the contributions from planetary and baroclinic waves, especially jet stream loops with typical spacings of 2000 km. Note that like higher altitudes, the levels of the structure function

_{LL}*D*and

_{LL}*D*merge at the smallest spacings. This behavior is consistent with the results of Cho and Lindborg (2001, their Table 2 and their Fig. 5) for the stratosphere.

_{NN}## 4. Results

In this section we present some results of the velocity and temperature structure functions (using the parameters in Tables A1 –A10 of the appendix) for various global and continental U.S. (CONUS) latitude bands and for various altitudes, seasons, and geographic regions. In Eq. (10), *s* has units of kilometers and, consequently, *a*_{1} has units of m^{2} s^{−2} km^{−2/3}, *a*_{2}, *a*_{4}, and *c _{i}* have units of kilometers,

*a*

_{3}is unitless, and

*b*has units of meters squared per second squared. Note that the empirical models for temperature structure functions in Tables A1 and A2 are produced from the AMDAR data only since the results from the sounding data have poor accuracy because of the lack of data in the tropics. We also derive vertical profiles of the turbulence levels,

_{i}*ε*

^{2/3}and

*C*

_{T}^{2}, and a geographic distribution of

*ε*

^{2/3}at upper levels (250 hPa) in regions where data coverage is adequate. Finally, as an independent check, the velocity and temperature structure functions are compared with the 20-km Rapid Update Cycle (RUC20) NWP-derived structure functions for a limited altitude and latitude band.

### a. Structure-function results versus altitude

Examples of the velocity and temperature structure functions at various pressure levels for the 40°–50°N latitude band are shown in Fig. 9. The largest variances are at the 200–300-hPa levels at nearly all scales corresponding to the jet stream maximum at roughly 200 hPa (Koch et al. 2006, their Fig. 1). The mesoscales have similar power but different slopes as a function of altitude. The velocity structure functions show a marked transition in slope at ∼100-km separation, consistent with Cho and Lindborg (2001). The smallest scales consistently show the *s*^{2/3} behavior at all altitudes. For lower altitudes the structure function levels decrease with decreasing altitude for the larger separations *s*. These results are consistent with Gage and Nastrom (1986). The larger scales show a consistent peak at all altitudes at separations in the 3000–10 000-km range for *D _{LL}*, with a narrower, more pronounced peak in

*D*at ∼2000 km or a wavelength of ∼4000 km. At this latitude band, the circumference of the earth is approximately 28 000 km, thus the broad maximum in

_{NN}*D*corresponds to planetary waves of zonal wavenumbers 3–5 while the peak in

_{LL}*D*corresponds roughly to zonal wavenumber 7 in agreement with the results of Benton and Kahn (1958).

_{NN}However, *D _{T}* shows a somewhat different behavior. At the largest lags the levels increase with decreasing altitude, except for a large increase at 200 hPa. At smaller lags there is no systematic change in level with altitude. From Hoinka (1998), the annual mean tropopause level for the 40°–50°N latitude band is ∼250 hPa. Thus

*D*at 200 hPa is mainly in the stratosphere and the increased structure-function level there is consistent with Nastrom and Gage (1985; Tables A1, A2). The transition point is at slightly larger scales (200–300 km) in

_{T}*D*compared with

_{T}*D*and

_{LL}*D*. There is no obvious peak in

_{NN}*D*at larger scales, just a slow roll-off starting at ∼1000–5000 km. At 500 hPa the roll-off starts at approximately 4000 km or wavenumber 7, in agreement with Julian and Cline (1974) at 40°N.

_{T}### b. Structure functions—Latitudinal variations

Figure 10 shows the latitude dependence of *D _{LL}*,

*D*, and

_{NN}*D*at 250 hPa. The tropical regions below 20°N have the lowest power at all scales, especially for temperature. The power generally increases with latitude, especially for the larger scales, although for spatial scales between approximately 100 and 1000 km the velocity structure functions are nearly independent of latitude above ∼30°N. The maximum of

_{T}*D*is in the 40°–50°N latitude range for all scales. The peak in

_{T}*D*is observed only in the northern latitudes above 30°N. All these results are consistent with the Nastrom and Gage (1985, their Fig. 8) GASP spectra and Cho and Lindborg (2001, their Fig. 2) velocity structure functions from the MOZAIC tropospheric data. Since our analysis provides information at larger scales, a zonal wavenumber 3 peak is also apparent in the 30°–40°N latitude range. This is consistent with three geographically forced jet maxima climatologically centered at approximately 80°W, 40°E, and 120°E in roughly the 30°–40°N latitude band (Koch et al. 2006, their Fig. 1). However, this feature is only noticeable in

_{NN}*D*. Since it is not evident in

_{NN}*D*, this apparent wavenumber 3 feature could be an artifact of the concentration of data over the continents.

_{LL}All latitude regions show the ubiquitous *s*^{2/3} scaling at small separations at this altitude (250 hPa). The transitional scale where the slope of the structure function increases from ⅔ is at ∼50 km, and is only weakly dependent on latitude except in the tropics where there is no obvious increase in slope. Nastrom and Gage (1985) and Cho and Lindborg (2001) also show only a weak dependence on latitude for the transition scale away from the tropics.

### c. Structure functions—Seasonal variations

The global structure functions for the 250-hPa pressure level for the summer (June–August) and winter (December–February) seasons for 40°–50°N and 30°–40°N are shown in Figs. 11 and 12. In general, the winter has higher levels at all scales, but this enhancement is especially noticeable in the velocity structure functions at 30°–40°N. The higher levels in this latitude band are consistent with the altitude structure at 40°–50°N described in section 4a and are related to the geographical distribution of jet streams, which is most pronounced in the winter months (Grotjahn 1993, his Fig. 5.7; Koch et al. 2006, their Fig. 4). Enhanced variability at the smaller scales is also observed in the winter months based on pilot reports of turbulence from aircraft (Wolff and Sharman 2008, their Fig. 5) and from global model reanalysis data (Jaeger and Sprenger 2007, their Figs. 1, 2).

Most noticeable is that the peak in *D _{NN}* at approximately 2000-km separation (zonal wavenumber 7) is apparent only in the winter, reinforcing the notion that this is related to planetary free waves (Benton and Kahn 1958). At both latitudes shown,

*D*also has an enhanced component at 8000 km (zonal wavenumber 3) in the winter, which is most likely related to the wavenumber 3 feature noted above.

_{LL}### d. Structure functions—Continental variations

The global velocity and temperature structure functions for 250 hPa and the 40°–50°N latitude band for the continental United States (50°–140°W) and Eurasia (10°W–160°E) are shown in Fig. 13. There is little difference in the temperature statistics at larger scales, but the smaller scales have substantially lower levels over Eurasia than over CONUS. The levels of *D _{LL}* and

*D*also have higher levels over CONUS than over Eurasia at the smaller scales, probably due to enhanced mountain wave activity over the Rocky Mountains in winter and to thunderstorm activity over the eastern half of the United States in summer (Wolff and Sharman 2008), compared to Eurasia. At larger scales

_{NN}*D*has the zonal wavenumber 7 peak mentioned in section 4c, and dominates

_{NN}*D*over both regions. At the largest scales, the CONUS levels are again higher than the Eurasia levels.

_{LL}### e. Structure functions—Effects of mountains

The impact of mountainous regions is investigated by selecting two subdomains with a high density of AMDAR data: the Rocky Mountains region (latitude 30°–50°N, longitude 100°–125°W) and the eastern United States (latitude 30°–50°N, longitude 75°–100°W). The comparison of the velocity structure functions is shown in Fig. 14. Note the large enhancement for the mountain (western) region at lags less than 20 km as well as the increase in power of *D _{LL}* compared to

*D*. This behavior has also been observed by similar analysis of research aircraft data (Frehlich and Sharman 2008), and is consistent with Jasperson et al. (1990) and the predictions of elevated turbulence spectra levels dominated by gravity waves (Cho et al. 1999; Lindborg 2007).

_{NN}### f. Profiles of small-scale turbulence

Recent theoretical and experimental results for stably stratified anisotropic turbulence (Lindborg 1999, 2006; Cho and Lindborg 2001; Riley and Lindborg 2008) have demonstrated that the *s*^{2/3} scaling of *D _{LL}*,

*D*, and

_{NN}*D*extends to the smallest scales where the theory of homogeneous isotropic turbulence is valid [Tatarski 1967, his Eqs. (2.7) and (3.19), respectively]. Then,

_{T}where *s* is in meters and *C* ≈ 2.0 is the Kolmogorov constant (Frisch 1995, p. 90). The AMDAR data resolve spatial scales of 10–50 km with sufficient statistical accuracy to produce a robust *s*^{2/3} behavior. Consequently, Eqs. (13) and (14) can be used to derive reliable estimates of the average energy dissipation rate *ε* and structure constant *C _{T}*

^{2}. Third-order structure functions can also be used to estimate

*ε*but they have much larger statistical estimation error.

Comparing Eq. (13) with Eq. (10) gives *ε*^{2/3} = *a*_{1}/100*C* = *a*_{1}/200 from *D _{LL}* and comparing Eq. (14) with Eq. (10) gives

*C*

_{T}^{2}=

*a*

_{1}/100 from

*D*. The coefficients

_{T}*a*

_{1}(m

^{2}s

^{−2}km

^{−2/3}) are listed in Tables A1 –A10 of the appendix for various global and CONUS latitude bands and for various pressure levels. As an example, the vertical profiles of meridionally averaged

*ε*

^{2/3}and

*C*

_{T}^{2}are shown in Fig. 15 for the 40°–50°N latitude band. A minimum in both the thermal (

*C*

_{T}^{2}) and velocity (

*ε*

^{2/3}) turbulence levels occurs at ∼8-km altitude, which reflects the lower average turbulence conditions above boundary layer convection and below the higher turbulence conditions associated with enhanced shear near the jet stream. Above this low-turbulence region the turbulence levels steadily increase to ∼200 hPa, near the ceiling of most commercial aircraft. Comparison of the

*a*

_{1}values in the tables of the appendix for other regions shows this vertical distribution is typical. Note that the structure functions at the lower altitudes (pressure greater than approximately 400 hPa; see Fig. 9) have a smaller

*s*

^{2/3}regime than the Lindborg model based on typical aircraft cruising altitudes of 9–11 km. The scatter in the data at lower altitudes in Fig. 15 are produced by the smaller number of observations at the small lags and from the contribution from random outliers, which are difficult to identify because of the large tails in the distribution of the velocity and temperature differences (Figs. 3, 4).

These values of *ε* and *C _{T}*

^{2}are consistent with previous global estimates [e.g., in Fig. 15 at 10-km elevation,

*ε*

^{2/3}≈ 1.4 × 10

^{−3}m

^{4/3}s

^{−2}or

*ε*≈ 5.2 × 10

^{−5}m

^{2}s

^{−3}, which agrees with

*ε*≈ 7.64 × 10

^{−5}m

^{2}s

^{−3}from Lindborg (1999) and

*ε*

^{2/3}≈ 1.7 × 10

^{−3}m

^{4/3}s

^{−2}for the 30°–50°N latitude region from Cho and Lindborg (2001, their Table 1)]. Similarly, at 10 km,

*C*

_{T}^{2}≈ 4 × 10

^{−4}K

^{2}m

^{−2/3}agrees with

*C*

_{T}^{2}≈ 5 × 10

^{−4}K

^{2}m

^{−2/3}determined by Frehlich and Sharman (2004) using the Nastrom and Gage (1985) GASP temperature spectra.

An early climatology of *ε* by Ellsaesser (1969) has a similar altitude dependence but the values are higher, most likely a consequence of using large lags around 200 km that are not in the *s*^{2/3} scaling region. Vinnichenko and Dutton (1969) present velocity spectra in the free atmosphere with the low-turbulence regime consistent with Fig. 15, but the high-turbulence values of *ε* are in excess of 0.01 m^{2} s^{−3} (corresponding to “severe CAT”). Similar variations in *ε* were reported by Chen (1974). Profiles of *ε* from ground-based Doppler lidar have free-atmosphere values around *ε* ≈ 1 × 10^{−5} m^{2} s^{−3} (Frehlich et al. 2006; Frehlich and Kelley 2008), which agree with Fig. 15, but boundary layer values of *ε* can vary from as low as 1 × 10^{−6} m^{2} s^{−3} for stable nighttime conditions to approximately 0.01 m^{2} s^{−3} for convection (Frehlich et al. 1998, 2006) to as much as 0.1 m^{2} s^{−3} for high wind conditions (Frehlich and Kelley 2008). Still another source for comparisons is the climatology of turbulence derived from NWP model output, which predicts 〈*ε*^{2/3}〉 = 1.7 × 10^{−3} m^{4/3} s^{−2} for an altitude of 10 km (Frehlich and Sharman 2004), which again is in good agreement with Fig. 15.

Other three-dimensional estimates of turbulence climatology have been attempted by Jaeger and Sprenger (2007) using 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) data and by Wolff and Sharman (2008) using pilot reports (PIREPS) of turbulence. However, both of these techniques are somewhat unreliable given the uncertainty in turbulence diagnostic formulations and model data used in the former and the subjective nature of PIREPs in the latter. Neither technique provides quantitative estimates of turbulence intensity.

Another point of comparison is inferences of *ε* and *C _{T}*

^{2}from radar measurements of spectral width or profiles of refractivity (e.g., Gage et al. 1980; Fukao et al. 1994; Cohn 1995; Hocking 1996; Nastrom and Eaton 1997, 2005; Satheesan and Krishna Murthy 2002; Masciadri and Egner 2006) or from stability profiles derived from high-resolution rawinsonde data (e.g., Bertin et al. 1997; Clayson and Kantha 2008). However, these are usually based on individual cases and tend to have large scatter.

Climatological studies of radar-derived *ε* and *C _{T}*

^{2}have been performed, but they are tied to a particular geographic location. Longer records of radar-derived

*ε*(4 yr in the United States) by Nastrom and Eaton (1997) over White Sands, New Mexico, and Nastrom and Eaton (2005) over Vandenberg Air Force Base, California, have similar shapes. Similar results were produced by Fukao et al. (1994) over Japan. In comparing these results, it should be noted that the inference of

*ε*from radar statistics is not straightforward, and involves a number of assumptions (Hocking 1996; Cohn 1995) regarding the nature of atmospheric turbulence in the free atmosphere, which may increase the uncertainty associated with radar-derived turbulence metrics. Our estimates of

*ε*and

*C*

_{T}^{2}are more directly tied to actual wind and temperature measurements and have three-dimensional coverage, and therefore should be more reliable and more indicative of true climatologies than those provided by limited radar studies.

From Eq. (10), the level of the structure functions at small spacings *s* scales as *a*_{1}*s*^{2/3}, where *a*_{1} represents *a _{L}*,

*a*, and

_{N}*C*

_{T}^{2}for the longitudinal velocity, transverse velocity, and temperature structure functions, respectively. An important aspect the observed levels in assessing theories of mesoscale turbulence is the ratio of the kinetic energy spectrum Φ

_{KE}(

*k*) to the potential energy spectrum Φ

_{PE}(

*k*) (Gage and Nastrom 1986; Gage et al. 1986; Cot 2001; Lindborg 2005, 2006; Lindborg and Brethouwer 2007). For the higher wavenumbers with the

*k*

^{−5/3}scaling corresponding to the

*s*

^{2/3}scaling of the structure functions (the spectral level is proportional to the structure function constant

*a*

_{1}) (Cot 2001),

is independent of *k*, where *T*_{0} is the average absolute temperature, *g* is the acceleration of gravity, and *N* ^{2} = (*g*/*θ*)(∂*θ*/∂*z*) where *θ* is the mean potential temperature, *z* is the vertical direction, and *N* is the Brunt–Väisälä frequency. Since vertical gradients cannot be estimated from the AMDAR data, *N*(*z*) is calculated from the U.S. Standard Atmosphere. The ratio *R* is tabulated in Table 1 for a few chosen altitudes in the troposphere and one in the stratosphere. Note that *R* is not constant with altitude but varies from approximately 1 to 2 with the highest values in the stratosphere (200 hPa). A similar enhancement in the stratosphere was produced by analyses of high-rate balloon data (Nastrom et al. 1997) but *R* was approximately a factor of 2 higher. Our magnitudes of *R* are consistent with other estimates from theory, observations, and numerical simulations of stratified and rotating turbulence. For example, *R* = 2 was proposed by Charney (1971) for geostrophic turbulence, and for plane inertia–gravity waves propagating in a constant wind and stability environment *R* ≥ 1, where the lower limit of unity corresponds to no rotation (Gill 1989; Nappo 2002). Gage and Nastrom (1986; Gage et al. 1986) found *R* ≈ 2 based on the GASP data. Lindborg (2006) produced a ratio (*R* = *ε _{K}*/

*ε*in his notation) that varied from 2.45 to 3.16 from numerical simulations of stratified turbulence. Including rotation in the simulations, Lindborg (2005) produced ratios that varied from 1.67 to 2.33. The ratios from Lindborg and Brethouwer (2007) varied from 1.7 to 3.1.

_{P}A map of average *ε*^{2/3} determined from all the AMDAR-derived velocity structure functions for 10° latitude–longitude boxes at 250 hPa for the Northern Hemisphere is shown in Fig. 16. There is a pronounced maximum over the Rocky Mountains (see Fig. 14), consistent with previous results (Nastrom and Gage 1985; Jasperson et al. 1990; Jaeger and Sprenger 2007; Wolff and Sharman 2008). The maximum off the southeast coast of the United States was also found in Wolff and Sharman (2008). The other elevated turbulence region over Asia is consistent with Jaeger and Sprenger (2007).

### g. Comparison with NWP model output

Average velocity and temperature structure functions derived from the archive of the RUC13 analyses (interpolated to the RUC20 grid) for the years 2005 and 2006 are compared with the soundings/AMDAR average structure functions calculated over the same spatial domain in Fig. 17. The agreement between the two estimates for scales larger than 400 km is better than 10%, which indicates a robust climatology with just 2 yr of NWP model output. Similar results were produced from a 1-yr analysis of the Global Forecast System (GFS) output (not shown).

## 5. Summary and discussion

Structure functions [Eq. (1)] are easily computed from routine measurements of wind and temperature provided by the global AMDAR and rawinsonde observations, and allow computation of turbulence levels without recourse to data provided through special sampling programs such as GASP and MOZAIC. The large archive of rawinsonde and AMDAR data samples the 5–10 000-km spatial scales of atmospheric processes with acceptable accuracy for northern latitudes. Spatial scales less than 100 km are the most challenging since small observation errors or random outliers can have a large impact on the long-term averages (see Figs. 3, 4) and therefore careful QC of the data is required. Nevertheless, structure function shapes and levels were found to be very robust for scales as small as 2 km for velocity and 10 km for temperature. It was found that the statistics from the sounding data are more accurate for the larger scales while the AMDAR data provide better estimates at the smaller scales. Good agreement is produced in the overlap scales of 200–500 km (Figs. 5, 6) even though the spatial and temporal averages of the two data sources are quite different. Thus the combination of rawinsonde/AMDAR data to derive structure functions produces a complete description for all scales greater than ∼5 km (see Figs. 7, 8). Empirical models [Eq. (10)] are fit to the merged structure functions for all regions that have sufficient statistical accuracy (see appendix).

Examples of the altitude dependence of the structure functions are shown in Fig. 9 and of latitude dependence in Fig. 10. The Lindborg model for the velocity structure functions is in good agreement with the results for the 250-hPa pressure level and the 40°–50°N latitude band (Fig. 7), which contains the largest number of commercial aircraft flights. There is a noticeable seasonal dependence at 30°–40°N and 250 hPa (Fig. 12), with more energy at the larger scales in the winter season. For the same region, there is good agreement at the mesoscales for the CONUS and Eurasian continent but with different behavior of the velocity statistics at the planetary scales. The different behavior of the structure functions at larger scales (greater than 2000 km) is partly related to the high density of data over the continents (Fig. 1) and therefore does not produce a true global average.

The leading-order scaling *s*^{2/3} of the empirical model for the structure functions of velocity and temperature can be used to produce profiles of the small-scale turbulence statistics *ε*^{2/3} and *C _{T}*

^{2}. The minimum of the turbulence levels at an altitude of 8 km (400 hPa) for the 40°–50°N latitude band reflects the lower turbulence region below the jet stream and above the convective boundary layer. This climatology of smaller-scale turbulence may provide useful connections to scales of turbulence that impact aviation safety (Wolff and Sharman 2008).

Our derived profiles of *ε*^{2/3} and *C _{T}*

^{2}predict the ratio

*R*of kinetic energy to potential energy in the

*k*

^{−5/3}region for the spatial spectra [Eq. (15); Gage and Nastrom 1986; Gage et al. 1986; Cot 2001; Lindborg 2005, 2006; Lindborg and Brethouwer 2007], assuming the Brunt–Väisälä frequency

*N*is a constant determined from the U.S. Standard Atmosphere. The ratio varied from 1.0 to 2.0, which is consistent with some past results and theoretical predictions, but the assumption of constant

*N*requires more investigation since local regions of high turbulence would be expected to produce small or even negative values of

*N*

^{2}(e.g., Bertin et al. 1997), which could bias the ratio

*R*. Including the effects of random variations in

*N*requires an appropriate spatial averaging domain for estimates of the

*θ*gradients (Reiter and Lester 1968; Balsley et al. 2007).

The climatology of turbulence will improve as more and better data are archived, especially for the smaller scales, which are most sensitive to random measurement errors and outliers. This climatology is essential for a rigorous evaluation of the effective spatial resolution of NWP models (Frehlich and Sharman 2004, 2008), for correct interpretation of forecast error statistics or innovation statistics (Frehlich 2008), and for optimal data assimilation (Frehlich 2006; Frehlich and Sharman 2004). Further, knowledge of the climatology of turbulence permits more accurate calculations of the total observation errors and forecast errors. Total observation errors consist of two components (Frehlich 2001): the instrument error and the observation sampling error (related to the error of representativeness), which defines the errors produced by the mismatch between the observation sampling volume and the definition of “truth.” The most precise definition of truth for a given NWP model is the convolution of the continuous atmospheric field by the spatial filter of the NWP model (Frehlich 2006, 2008). Rawinsonde data have the largest observation sampling error since the sampling volume is very small (i.e., it is approximately a single point in the 2D horizontal plane). The average structure function and the spatial filter of the NWP model then determine the average observation sampling error (Frehlich 2001) as well as the atmospheric contribution to the average forecast errors, also called innovation errors (Frehlich 2008). For example, considering the latitude region 40°–50°N and a global model with grid resolution Δ = 35 km and effective model resolution *L* = 150 km, the observation sampling error at 250 hPa for rawinsonde observations at the center of a grid cell for one horizontal velocity component is 1.61 m s^{−1} [Frehlich 2001, his Eq. (93)], which is larger than the instrument error of ≈0.5 m s^{−1}. In contrast, the observation sampling error for a rawinsonde temperature measurement is 0.50 K, which is comparable to the instrument error of 0.5 K.

Last, note that, as more and more aircraft become AMDAR equipped, it may be feasible to compute structure functions locally in real time and thus infer locations of elevated turbulence levels and provide new tools for turbulence detection and aviation safety (Wolff and Sharman 2008).

## Acknowledgments

We thank Bill Moninger, NOAA/ESRL/GSD, for supplying us with the quality controlled AMDAR data and for his helpful comments on the data. We also thank Teddie Keller for useful comments and suggestions on an earlier version of the manuscript. This work was funded in part by NSF Grants ATM0646401 and ATM0522004, ARO Grant W911NF-06-1-0256, and by NASA ROSES Grant NNX08AL89G. We appreciate the very helpful comments of Erik Lindborg and the anonymous reviewers.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

^{3}and k

^{5/3}energy spectrum of atmospheric turbulence: Quasigeostrophic two-level model simulation.

**,**

**,**

**,**

**,**

### APPENDIX

#### Coefficients of Best-Fit Models

The coefficients of the best-fit models [Eq. (10)] for the global-averaged structure functions are given in Tables A1 –A10 for the Northern Hemisphere using the merging technique shown in Figs. 7 and 8, except for the temperature structure functions in the latitude bands less than 20°N, which are based on the AMDAR data only. In Eq. (10), the variable *s* has units of kilometers, and consequently *a*_{1} has units of m^{2} s^{−2} km^{−2/3}, *a*_{2}, *a*_{4}, and *c _{i}* have units of kilometers,

*a*

_{3}is unitless, and

*b*has units of meters squared per second squared.

_{i}## Footnotes

*Corresponding author address:* Rod Frehlich, CIRES, UCB 216, University of Colorado, Boulder, CO 80309. Email: rgf@cires.colorado.edu

* The National Center for Atmospheric Research is sponsored by the National Science Foundation.