ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0 FROM SURFACE-LAYER WIND DATA D. ANFOSSI1 , G. DEGRAZIA2 , E. FERRERO3 , S. E. GRYNING4 , M. G. MORSELLI5 and S. TRINI CASTELLI1 1 CNR-Istituto di Cosmogeofisica, Torino, Italy; 2 Universidade Federal de Santa Maria (RS), Brasil; 3 Universita’ Piemonte Orientale, Dip. Scienze Tecn. Avanz., Alessandria, Italy; 4 RISO National Laboratory, Roskilde, Denmark; 5 ENEL/SRI Area Ambiente, Milano, Italy

(Received in final form 27 December 1999)

Abstract. Eulerian turbulence observations, made in the surface layer under unstable conditions (z/L < 0), by a sonic anemometer were used to estimate the Lagrangian structure function constant C0 . Two methods were considered. The first one makes use of a relationship, widely used in the Lagrangian stochastic dispersion models, relating C0 to the turbulent kinetic energy dissipation rate , wind velocity variance and Lagrangian decorrelation time. The second one employs a novel equation, connecting C0 to the constant of the second-order Eulerian structure function. Before estimating C0 , the measurements were processed in order to discard non-stationary cases at least to a first approximation and cases in which local isotropy could not be assumed. The dissipation  was estimated either from the best fit of the energy spectrum in the inertial subrange or from the best fit of the third-order longitudinal Eulerian structure function. The first method was preferred and applied to the subsequent part of the analysis. Both methods predict the partitioning of C0 in different spatial components as a consequence of the directional dependence of the Eulerian correlation functions due to the isotropy. The C0 values computed by both methods are presented and discussed. In conclusion, both methods provide realistic estimates of C0 that compare well with previous estimations reported in the literature, even if a preference is to be attributed to the second method. Keywords: Lagrangian and Eulerian statistics, Structure functions, Autocorrelation functions, Energy spectra, Surface layer.

1. Introduction The atmospheric dispersion of airborne pollutants can be successfully simulated with Lagrangian stochastic dispersion models, based on the generalised Langevin equation (Wilson et al., 1983; Thomson, 1987; Sawford and Guest, 1988; Luhar and Britter, 1989; Weil, 1990; Flesh and Wilson, 1992; Wilson and Flesh, 1993; Du et al., 1994; Rotach et al., 1996; Wilson and Sawford, 1996; Ferrero and Anfossi, 1998a,b). In these models (Ito type models), it is assumed that the velocity ui (i = 1, 2, 3: along-wind, across-wind, and vertical) of a fluid element, for times larger than the Kolmogorov time scale, is a Markov process that is a continuous function of time (Sawford, 1993). The equations prescribing the evolution of ui and of the displacement xi are (Thomson, 1987; Wilson and Sawford, 1996) p dui = ai (x, u, t) dt + C0 δij dξi (t), (1a) Boundary-Layer Meteorology 95: 249–270, 2000. © 2000 Kluwer Academic Publishers. Printed in the Netherlands.

250

D. ANFOSSI ET AL.

and dxi = ui dt,

(1b)

where ai (x, u, t) is a coefficient that is determined from the Fokker–Planck equation for stationary condition and depends upon the Eulerian probability density function P (x, u), C0 is a numerical constant,  is the ensemble-average rate of dissipation of turbulent kinetic energy and dξi (t) is a component of Gaussian white noise, which is uncorrelated with other components and is uncorrelated in time, i.e., hdξj (t) dξi (t + s)i = δij δ(s) dt ds.

(2a)

The two terms on the r.h.s. of Equation (1a) represent the deterministic and random velocity forcing, respectively. To derive the expression of the random forcing √ term, C0 δij (Van Dop et al., 1985; Thomson, 1987; Rodean, 1994; Wilson and Sawford, 1996), let us write Equation (1a) in its more general form dui = ai (x, u, t) dt + bij (x, u, t) dξi (t),

(2b)

where bij (x, u, t) is a coefficient to be determined. Then, the turbulence structure function (the ensemble average of the square of the change in Lagrangian velocity) can be evaluated directly for small times from Equations (2a and b) and written here as DijL (τ ) = bij2 τ, where τ is an infinitesimal time increment in the Kolmogorov inertial subrange. By equating this to the second-order Lagrangian velocity structure function, determined according to Kolmogorov’s theory of locally isotropic turbulence in the inertial subrange (Monin and Yaglom, 1975), DijL (τ ) = C0 δij τ, we obtain p bij = C0 δij . In particular, this means that the Lagrangian structure function constant C0 is an important quantity in model simulations. However, the value of C0 is not yet accurately established and still represents a matter of discussion (Rodean, 1994). Hinze (1975) found that C0 varies between 3 and 10. A fit of modelled vertical dispersion parameters σz from a line source in the constant stress region of a neutral boundary layer indicates that the value of C0 lies between 5.0 and 10.0 (Sawford and Guest, 1988). Hanna (1981) experimentally evaluated C0 = 4.0 ± 2.0 in the convective boundary layer (CBL), Anand and Pope (1985) found good results in

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

251

wind-tunnel measurements by using C0 = 2.1, whereas Sawford and Tivendale (1992) and Sawford and Borgas (1994) found that C0 = 3 gave the best fit. Du et al. (1995), by comparing measured and modelled dispersion in decaying turbulence in a water channel and wind tunnel, found 3.0 ± 0.5. Horst and Weil (1992) used C0 = 5 to estimate footprints for scalar fluxes in the convective surface layer. Theoretical considerations yield C0 = 7 (Sawford, 1991) and C0 = 5.7 (Rodean, 1991). Degrazia and Anfossi (1998) theoretically derived C0 = 3.4 for the crosswind and vertical turbulent velocity components and C0 = 2.2 for the longitudinal one, in the CBL. For dispersion studies with stochastic models, a relative small value of C0 (between 2 and 3) is used by Luhar and Britter (1989), Hurley and Physick (1991, 1993), Physick et al. (1994), Tassone et al. (1994) and Rotach et al. (1996). This large variability in C0 estimates may imply that C0 is not a universal constant, supporting the suggestion of Landau and Lifshitz (1989) that C0 is a function of the type of flow. This argument was also treated at length by Frisch (1995). The original Landau and Lifshitz (1989) argument is here reported, only changing the notation according to the present paper and considering time scales instead of length scales. ‘It might be thought that the possibility exists in principle of obtaining a universal formula applicable to any turbulent flow which should give Dil (τ ) for all times τ (that are small compared with TiL ). The instantaneous values of DiL (τ ) might in principle be expressed as a universal function of the dissipation  at the instant considered. When we average these expressions, however, an important part will be played by the manner of variation of  over times of the order of the periods of the large eddies, and this variation is different for different flows. The result of the averaging therefore cannot be universal’. We considered C0 as a nonuniversal constant, following the Landau and Lifshitz (1989) and Frisch (1995) suggestions. Consequently, C0 can be estimated in terms of other non-universal constants. In this paper, two alternative experimental evaluations of C0 , based on wind data recorded by a sonic anemometer in the atmospheric surface layer under unstable conditions, are presented and discussed. In Section 2 the theory used in this work is briefly presented, Section 3 describes the site and the data, Section 4 deals with the data analysis, and discussion of results appears in Section 5. 2. Theory Two methods for the determination of C0 from measured wind velocity fluctuations are proposed. In the first method (M1), the following relationship is used C0 =

2σi2 TiL

t < TiL ,

(3)

where σi2 is the velocity fluctuation variance and TiL is the Lagrangian decorrelation time scale. Equation (3), which is widely used in Lagrangian stochastic

252

D. ANFOSSI ET AL.

models, was derived by Hinze (1975) and Tennekes (1982), assuming an exponential form for the autocorrelation function. The calculation of σi2 from measured wind velocity fluctuation data is straightforward. The eddy dissipation rate  may be computed from the inertial subrange of the observed Eulerian wind spectra (Hanna, 1981) by a best fit technique applied to the observed spectrum and the theoretical one SiE (n) = Ai (u) ¯ 2/3  2/3n−5/3 ,

(4)

in the inertial subrange (n being the frequency). The obtained  value depends on both the values of the empirical constants Ai and on the accuracy of the spectrum and its slope estimation. Another method is based on the third-order longitudinal Eulerian velocity structure function (Frisch, 1995) 4 DE3 (τ ) = h[u(t + τ ) − u(t)]3 i = −  uτ, ¯ 5

(5)

for increments τ small compared to the integral scale. In Equation (5), u¯ is the mean wind speed. In this case there are no empirical constants, however the determination of  may be less accurate due to the highly noisy character of the third-order structure function estimated from field measurements. The procedure to estimate the Lagrangian time scale TiL is as follows. The Eulerian autocorrelation function Ri (τ ) is estimated first. This may be done in various ways – from the best fit of Ri (τ ) assumed exponential, from the integration of Ri (τ ), and from the estimation of the time lag at which Ri (τ ) firsts falls to 0.37 (1/e). Finally, TiL is derived from the ratio between the Lagrangian and Eulerian integral time scales (Gifford, 1955; Hay and Pasquill, 1959; Wandel and Kofoed (1962); Corrsin, 1963; Angell et al. (1971); Csanady, 1973; Hinze, 1975; Hanna, 1981; Panofsky and Dutton, 1984; Sorbjan, 1989), TiL /TiE = βi = γ u/σ ¯ i,

(6)

where γ is a numerical coefficient. Though a complete theoretical derivation of Equation (6) is not yet available, such a relationship suggested by Gifford (1955) and Hay and Pasquill (1959) to interpret their experimental results was confirmed by other works (Wandel and Kofoed, 1962; Corrsin, 1963; Saffman, 1963; Smith, 1967; Slade, 1968; Angell et al., 1971; Hanna, 1981; Panofsky and Dutton, 1984) and thus we believe it can be used with a sufficient level of confidence. Using Equation (6), Equation (3) can be rewritten as C0i =

2σi3 1 ,  uT ¯ iE γ

(7)

in which C0i has been indicated instead of C0 in order to allow for different values according to the different wind components.

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

253

The second method (M2) is based on the calculation of the second-order Eulerian velocity structure function DiE (τ ); the procedure is as follows (see also Appendix A). Recall (Monin and Yaglom, 1975) that the expressions for the Lagrangian velocity structure function DiL (τ ) and the Lagrangian energy spectrum E L (ω), where ω is the angular frequency, in the inertial subrange, are the following, DiL (τ ) = C0i τ

(8a)

EiL (ω) = B0i ω−2 ,

(8b)

where B0i =

C0i π

(8c)

is a numerical constant. The corresponding Eulerian functions are expressed as DiE (τ ) = αi Cs ( u) ¯ 2/3 τ 2/3

(9a)

EiE (ω) = αi αu ( u) ¯ 2/3ω−5/3 ,

(9b)

where Cs αu ∼ . = 4

(9c)

Cs and αu are numerical constants and αi = 1, 4/3, 4/3 for u, v and w components, respectively, as a consequence of isotropy. Changing from angular frequency ω to frequency n, according to the Taylor’s ‘frozen flow’ approximation, i.e., setting ω = 2π n and S(n) = 2π E(ω), Equations (8b) and (9b) become SiL (n) =

B0i −2 n , 2π

(10a)

αi αu ( u) ¯ 2/3 n−5/3 . (2π )2/3

(10b)

and SiE (n) =

Pasquill (1974) showed that the ratio between the Lagrangian and Eulerian time scales can be expressed, with the present notation, as  3/2 TiL 3 (αi αu )3/2 u¯ = π . (11) E 2 B0i σi Ti Now, considering Equations (6), (8c) and (9c), we finally obtain   (Cs )3/2 γ 8 3/2 = , C0i π 3αi

(12)

254

D. ANFOSSI ET AL.

or C0u =

(Cs )3/2 , 1.39γ

(13a)

and C0v = C0w =

(Cs )3/2 . 0.90γ

(13b)

Cs is computed from the inertial subrange of the observed Eulerian structure function and the procedure to estimate  is the same as in M1. In conclusion, two methods are proposed to compute C0 from observed wind data (M1 and M2). For clarity reasons, C0 values obtained from M1 (Equation (7)) 1 will be named hereafter C0i and those computed from M2 (Equation (13)) will be 2 named C0i . It may be worth noting the similarities and the differences between the two approaches; in both there is the need of estimating  and of having a value for the coefficient γ to calculate C0i . M1 needs to compute the Eulerian autocorrelation function Ri (τ ) (for estimating the Eulerian time scale) and M2 needs to compute the structure function DiE (τ ) (to estimate the constant Cs ). However these computations are equivalent since the two functions are related by the following expression DiE (τ ) = 2[Ri (0) − Ri (τ )].

(14)

Therefore, the main difference is that M1 needs to estimate the Eulerian time scale TiE . Due to a certain degree of arbitrariness inherent in this determination, it might be expected that, in principle, M2 should be the more robust way of calculating C0i . Different component values of C0 (as in Equations (7) and (12)) were previously introduced by Degrazia and Anfossi (1998). It is appropriate to consider where this C0 directional partitioning originates from. In that paper, Degrazia and Anfossi (1998) starting from the present Equation (3), and accounting for the observed spectral properties and the classical statistical diffusion theory and making use of the above relationship between Lagrangian and Eulerian time scales (Equation (6)), obtained a theoretical expression that relates C0 to γ and αi . In the present M2 method, C0 is found to be a function of Cs , γ and αi . According to Monin and Yaglom (1975), the coefficient Cs (see Equation (9a)) should be constant. As a consequence, in M2 and in the Degrazia and Anfossi (1998) derivation, different C0 values for the different components derive from αi , that is from the fact that the Eulerian correlation functions are directionally dependent. In the M1 method, the C0 directional dependence derives from the ratio between σi2 and TLi (Equation (3)) or between σi2 and TEi (Equation (7)). Thus, in all these derivations, the C0

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

255

component partitioning originates from the chosen relationship between TLi and TEi . As mentioned above, this last relationship is based on the interpretation of experimental results and is widely quoted in the literature.

3. Wind Data Data used in the present analysis were gathered from the Northern Hemisphere Climate-Processes Land-Surface Experiment (NOPEX). The measurements were collected within agricultural fields of barley and wheat at the Tisby site situated about 45 km west of Uppsala, Sweden, at 16◦ 570 E and 59◦ 460 N. The area around Tisby is almost flat, though weakly sloping towards a small river. The typical distance to forested areas is 1.5 km in all directions except south, where the fetch is homogeneous for more than 10 km. The measurements were carried out in 1995 during a three-month period starting mid-April, and thus represent spring and summer conditions. The full NOPEX instrumentation and measurement aims and strategy are described in Thomsen et al. (1995). In this study we consider wind observations taken at a mast (6.8 m high) using an ultrasonic anemometer (type Kaijo Denki DAT/TR-61B) with a sampling frequency of 20 Hz.

4. Data Processing 27 days of measurements (from April 22 to May 17) were considered in the present analysis. The full data record was broken into 1-hour segments. Thus hourly raw data of wind (u, v, w) and temperature (T ) observations, each containing 72,000 data, were analysed. Hourly data sets are generally considered (Kaimal and Finnigan, 1994; Moore et al., 1985) long enough to obtain adequate estimates of the statistical quantities (autocorrelation, structure functions and spectra) computed in this work. The various steps of the data treatment, selection and computation are now presented. Wind data were firstly rotated to refer them to a streamline system (McMillen, 1988; Cassardo et al., 1995) and linear trends were removed from each series by the least squares method (Bendat and Piersol, 1971). Then the stationarity test suggested by Bendat and Piersol (1971) was applied to select stationary time series. The analysis only continued for those cases in which the three time series – u(t), v(t), w(t) – simultaneously passed the test. In the second phase, the wind variances σi2 and the surface-layer parameters (friction velocity u∗ temperature flux w 0 T 0 , Obukhov length L and scaling temperature T∗ ) were computed. Only those cases in which the mean horizontal wind speed was greater than 2 m s−1 were retained (to exclude the low wind speed cases in which the equations reported in Section 2 are not likely to be valid). Then,

256

D. ANFOSSI ET AL.

Eulerian structure functions DiE (τ ) and autocorrelation functions Ri (τ ) were computed by way of direct computation and velocity spectra SiE (n) were estimated by the standard method, i.e., by Ri (τ ) Fourier transform (Bendat and Piersol, 1971). In the third phase, Eulerian time scales TiE were estimated.Three alternative methods for computing TiE were introduced in Section 2: best fit of autocorrelogram, integration of autocorrelogram and determination of the time at which the autocorrelogram decreases to 1/e. In the first and third methods, an exponential form for Ri (τ ) is assumed. Hanna (1981) discussed the second and third method only. Since the second method was found to be less reliable, because the observed autocorrelogram does not generally approach zero at the largest time lags, Hanna chose the 1/e method. Moore et al. (1985) tested all three methods and found that the first one (best fit) was the more accurate. For this reason, in the present work we evaluated TiE from the best fit of Ri (τ ), assumed exponential. In the fourth phase, the ensemble-average rate of dissipation of turbulent kinetic energy  was estimated; two alternative methods for computing this quantity were considered (see Section 2). The first one consists in performing the best fit between the three wind velocity theoretical spectra (Equation (4)) and the corresponding observed ones, and the second one is based on a best fit of the third-order Eulerian longitudinal velocity structure function (Equation (5)). In both cases, only the inertial subrange is considered. In the first method, three values of  (u , v and w ) were computed from the three wind component spectra. That is, a best fit between the theoretical spectra (Equation (4)), in which Au = 0.15 and Av,w = 0.2 (Pasquill, 1974), and experimental spectra were performed, step by step, in the frequency range nmax < n < 1/TiE . Since the assumption of local isotropy (used in the derivation of Equation (4)) is likely to be only approximately satisfied in the real atmospheric surface layer, and due to a certain amount of arbitrariness in correctly locating the inertial subrange in the observed spectra and in estimating their exact slope, we cannot expect the three wind spectra to provide exactly the same value of . This point was also discussed by Kaiser and Fedorovich (1998) who found that, in general, their wind tunnel  values, computed for the longitudinal and vertical components, were rather different (see their Figure 9a). Consequently, to exclude the unrealistic cases, in this work a mean  value (¯ ) was computed, i.e., ¯ =

(u + v + w ) , 3

and only those cases in which all i satisfied the following relationships rel =

|i − ¯ | < 0.30 ¯

(15)

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

257

were selected. According to the Kolmogorov theory under conditions of local isotropy the spectral ratios Swu , Svu and Swv (where Sij = Si /Sj ), in the inertial subrange must be equal to 4/3, 4/3 and 1, respectively. To fulfil this basic requirement, at least approximately, a final selection was performed. Similarly to that done for , a mean value for each of the above ratios was computed and only those time series in which the following inequalities Srel =

|Swu − 4/3| |Svu − 4/3| |Svu − 1 < 0.30, = < 0.30, < 0.30, 4/3 4/3 1

were simultaneously satisfied, were considered for the final analysis. At the end of all the computation and selection procedures described so far, 77 1 2 cases (hours) remained. In the last phase, C0i and C0i values were estimated for these 77 selected cases, according to method M1 (Equation (7)) and method M2 (Equation (13)), by choosing γ = 0.55 according to Degrazia and Anfossi (1998). 5. Results Before presenting and discussing the comparison between the M1 and M2 estimates of C0i , it is necessary to choose between the two alternative ways of estimating  introduced in Section 2 and present some general results that will illustrate the overall quality of the estimates. To clarify which one of the two methods to determine  is the more reliable, the theoretically expected value (c ), according to Kaimal and Finnigan (1994), was computed for each time series from, c =

 z 2/33/2 u3∗ , 1.0 + 0.5 κz L

(16)

where κ is the von Karman constant. Figure 1 shows the comparison between these values (Equation (16)) and those estimated from the best fit of spectra in the inertial subrange; overall, there is fair agreement. The second approach (best fit of third-order longitudinal structure function, Figure 2) appears to be less reliable. As anticipated in Section 2, this result could be expected because the observed DE3 (τ ) are not very regular. This can be seen in Figure 3, where an example of observed −DE3 (τ ) versus τ is plotted (crosses). The straight line represents the −DE3 (τ ) ∝ τ 1 line, and refers to the Julian day 148 (28 May 1995 at 1000 hours). It appears that it is not easy to safely identify the region in which −DE3 (τ ) is linear in time. As a consequence, the ensemble average rates of dissipation of turbulent kinetic energy, calculated by means of the first approach (spectra), will be used in the remaining part of this work. The trends of DiE (τ ) versus time, SiE (n) versus frequency and Ri (τ )/Ri (0) versus time, are shown in Figures 4–6. They refer to the same case. In Figures 4 and 5 the straight lines represent the 2/3 and −5/3

258

D. ANFOSSI ET AL.

Figure 1. Comparison between turbulent kinetic energy dissipation rate  estimated from the spectral method (x-axis) and from Kaimal and Finnigan (1994) formula, Equation (16) (y-axis).

power law, respectively. It can be seen that the inertial subrange regions can be more easily identified in these two plots. The general quality of the plots obtained for the other cases is similar to the one represented by the chosen case. A test of the quality of the input wind data and of the data processing, at least in the first two phases, is illustrated in Figure 7. This plot shows the comparison between the present values of the standard deviation of the vertical component of the wind velocity, σw , as a function of z/L and those calculated from the following relationships (Garratt, 1992) σw = 1.25u∗ (1 − 3z/L),

(17)

which is widely accepted in the literature and used in practical applications (see, for instance, Zannetti, 1990). We did not compare the computed values of the horizontal wind variance components since, for these last quantities, Monin–Obukhov scaling does not apply (Garratt, 1992). An inspection of Figure 7 shows that the agreement is good. It is to be noticed that the numerical coefficient found in the present data is 1.20, i.e., very similar to 1.25 of Equation (17). The results depicted in Figures 1 and 7 suggest that the overall procedure is reliable. Since Equations (7) and (12) indicate that the estimate of C0i depend on 1/γ in both M1 and M2 methods, it is useful to rewrite those equations as follows, 1 C0i =

Fi1 , γ

with

Fi1 =

2σi3 ,  uT ¯ iE

(18)

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

259

Figure 2. Comparison between turbulent kinetic energy dissipation rate  estimated from the third-order longitudinal structure function, Equation (5), (x-axis) and from Kaimal and Finnigan (1994) formula, Equation (16) (y-axis). TABLE I 1 and F 2 = γ C 2 , where γ is the numerical constant Results of computation of Fi1 = γ C0i i 0i

appearing in the ratio between the Langrangian and Eulerian integral time scales (Equation 2 and C 2 are the Langrangian structure function constants, calculated by the (6)) and C0i 0i two methods M1 (Equation (18)) and M2 (Equation (19)). They represent mean values and their spreads over the 77 cases (each lasting 1 h) selected. These cases were recorded under unstable conditions (z/L < 0) in the atmospheric surface layer. Fu1

Fv1

Fw1

Fu2

Fv2

Fw2

1.78 ± 0.10

1.78 ± 0.11

1.99 ± 0.06

1.73 ± 0.02

2.36 ± 0.05

2.37 ± 0.02

F2 = i , γ

Fi2

and 2 C0i

 with

=π

3 αi C s 8

3/2 ,

(19)

and to present the results relative to Fi1 and Fi2 , which do not depend on any particular choice of γ . These are presented in Table I. It may be appropriate to first consider that the present estimate of the Cs value (see Equations (19) and (9a)), obtained for the three wind components, gave 1.8,

260

D. ANFOSSI ET AL.

Figure 3. Third-order Eulerian longitudinal structure function as a function of time lag τ (crosses). 3 (τ ) ∝ τ 1 line. The time series considered was recorded on 28 The straight line represents the −DE May 1995 at 1000 hours.

1.6 and 1.7 respectively, that is approximately constant (with an average value of 1.7), as implied in the Monin and Yaglom (1975) derivation. This is another consideration suggesting that the present analysis results may be reliable. Then, let us recall that the results shown in Table I were obtained from a sample of 77 different hourly time series, selected by means of the criteria discussed above, out of a period of 27 days. This suggests that the average values listed in Table I (and in the following Table II) are meaningful from the statistical point of view. Table I indicates that the agreement between Fu1 and Fu2 values is good, while the crosswind and vertical components exhibit a relative difference of 25% and 16%, respectively. It can be seen that the spreads of Fi1 are significantly greater than those of Fi2 . This probably reflects the fact that M1 has to estimate, beside the dissipation rate, the Eulerian time scale too. As mentioned above, and also discussed by Hanna (1981), this estimation shows a certain degree of arbitrariness and may be responsible for the greater spread. Furthermore, it may be noticed that the M2 estimate does not depend on the assumption of any functional form for the autocorrelation function, whereas M1 explicitly does. Finally, in order to calculate C0i we need a good and reliable estimate of γ . Degrazia and Anfossi (1998), considering a large number of theoretical and experimental works, computed an average value of 0.55. This value is also used here.

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

261

Figure 4. As in Figure 3, but for the second-order Eulerian longitudinal structure function, dotted-dashed line, dotted line and dashed line indicate u, v and w components, respectively. The straight line represents 2/3 power law.

TABLE II Results of computation of the Langrangian structure function constants C0 according to the two methods proposed in this work: C0i refer to M1 (Equation (7)) 2 to M2 (Equation (13)). The results obtained and C0i from the suggested γ value are shown in the first row. For comparison, C0 resulting from two extreme γ values (from Table I of Degrazia and Anfossi, 1998) are also indicated. γ

1 C0u

1 C0v

1 C0w

2 C0u

2 C0v

2 C0w

0.55 0.80 0.40

3.2 2.2 4.5

3.2 2.2 4.5

3.6 2.5 5.0

3.2 2.2 4.3

4.3 3.0 5.9

4.3 3.0 5.9

262

D. ANFOSSI ET AL.

Figure 5. As in Figure 3, but for the Eulerian energy spectra versus frequency. Continuous line, dotted line and dashed lines indicate u, v and w components, respectively. The straight line represents −5/3 power law.

Figure 6. As in Figure 3, but for the Eulerian autocorrelation functions. Continuous line, dotted line and dashed lines indicate u, v and w components, respectively.

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

263

Figure 7. Comparison between σw computed from the sonic anemometer (x-axis) and computed from Equation (17) (y-axis).

The results are shown in Table II, first row. The reason for the different C0 values proposed in the literature may become clearer when we notice that the coefficient γ is not a universal constant but depends on the type of turbulent flow (Hinze, 1975). To account for this fact, and for comparison purposes, the second and third rows of Table II report the C0i values obtained by using the extreme γ values in Table I of Degrazia and Anfossi (1998): γ = 0.80 (Saffman, 1963) and γ = 0.40 (Angell et al., 1971). Instead of choosing a value for γ by averaging existing estimates, and using this to estimate C0 , an alternative approach would be to take experimental (literature) values of C0 itself and average. This would not be a better choice because the variability of C0 estimates (2–10) is much larger than that for γ (0.4–0.8). Referring to Table II, and taking account of Table I, suggests the following considerations: – the C0i values found with both methods, though including their possible range of variation indicated by their spread and the different γ values, fit well with the range of values suggested in the literature as summarised in the Introduction; – in particular, the different values proposed for γ may partly explain the different C0 values found in the literature; – while most of the previous C0 estimates assumed a single value for the three wind components, Degrazia and Anfossi (1998) theoretically predicted the

264

D. ANFOSSI ET AL.

same C0 value for the crosswind and vertical components and a lower value for the longitudinal component. This is a consequence of the different values of the constants in the Eulerian spectrum (Equation (4)). Present M2 results seem to support this prediction. On the contrary, M1, though confirming the prediction of different C0 values for different components, yields the same value for the two horizontal components; 1 1 – C0v and C0w values (first row) are in fair agreement with those theoretically 1 obtained by Degrazia et al. (1998), namely 3.44. The fact that C0i are more 2 similar to the values predicted by Degrazia et al. (1998) than C0i may not be surprising, since both estimations, the theoretical and the experimental ones, 1 are based on the same relationship (Equation (7)). On the contrary C0u and the corresponding predicted value (2.23) are different; – M2 values are in good agreement with the one (C0 = 4) experimentally obtained by Hanna (1981). It may be concluded that both methods, M1 and M2, are reliable but M2 is preferred.

6. Conclusions In this work an experimental evaluation of the Lagrangian structure function constant C0 , based on surface-layer wind data gathered by a sonic anemometer in a flat region in southern Sweden, is presented. At present, the exact value of C0 is still an open matter. Obtaining a correct C0 value, besides being important in fundamental atmospheric turbulence studies, it is also of great importance in the Ito type models for air pollution dispersion modelling. Two alternative ways of estimating C0 were proposed. The first, M1, is based on an equation established by Hinze (1975) and Tennekes (1982), (Equation (3)), relating the value of C0 to the turbulent kinetic energy dissipation rate , the wind velocity fluctuation variance σi2 and the Lagrangian decorrelation time scale TiL . Considering that, according to the Hay and Pasquill (1959) and Gifford (1955) hypothesis, the ratio between the Lagrangian and Eulerian integral time scales may be expressed as a function of the standard deviation σi , the mean wind speed u¯ and a numerical constant γ , the resulting Equation (7) is obtained. The second one, M2, is based on a novel relationship (Equations (12) or (13)) in which C0 depends upon the value of the numerical constant of the Eulerian second-order structure function DiE (τ ) (Equation (9a)), where τ is the time lag. In turn, this last depends on the turbulent kinetic energy dissipation rate  and the mean wind speed u. ¯ Both methods predict the partitioning of C0 in different spatial components, as previously derived by Degrazia and Anfossi (1998).This result is the consequence of using the relationship between TLi and TEi (Equation (6)) suggested by Hay and Pasquill (1959) and Gifford (1955) and confirmed or used by many other authors (see, for instance, Wandel and Kofoed, 1962; Corrsin, 1963; Saffman, 1963; Smith,

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

265

1967; Slade, 1968; Angell et al., 1971; Csanady, 1973; Hinze, 1975; Hanna, 1981; Panofsky and Dutton, 1984; Sorbjan, 1989) and of the directional dependence of the Eulerian correlation functions due to the isotropy. 27 days of continuous measurement, divided into hourly time series, were analysed. Wind data were firstly rotated to a streamline system and surface-layer parameters was computed. The analysis was focused on unstable cases (z/L < 0), and selected according to a stationarity criterion. Velocity variances and Eulerian autocorrelation and structure functions were computed, Eulerian time scales were estimated and from these the Lagrangian scales were determined (by means of Equation (6)). Then the velocity variance spectra were estimated from which, selecting only the parts belonging to the inertial subrange, the value of  was derived. An alternative method to estimate  was considered, which makes use of the thirdorder longitudinal Eulerian structure function. However, it became clear that this latter method gives less reliable results and, consequently, the spectral method was used throughout our study. Besides the stationarity criterion, three other selection criteria were used: cases in which u¯ was less than 2 m s−1 , or rel > 0.30, or Srel > 0.30 were excluded; 77 cases were thus selected. Finally, C0 was computed 1 2 with both methods M1 and M2, named C0i and C0i , and were compared, by using γ = 0.55 suggested by Degrazia and Anfossi (1998). It was found that our present C0 values, from both M1 and M2, are in good 1 agreement with current theoretical and experimental estimations. In particular C0i are more similar to those predicted by a recent theoretical derivation of C0 (De2 grazia and Anfossi, 1988), whereas C0i are in good agreement with a previous experimental estimation of Hanna (1981). Though all the main parameters, , Cs and T E , were estimated by means of a best-fit technique, thus introducing some arbitrariness into our results, it appeared that the third parameter is the most critical. This is reflected in the M2 results. Due to the different values of the constants in the Eulerian spectra (because of the local isotropy existence in the inertial subrange) a lower C0 value is expected 2 1 for the longitudinal wind component. While C0i confirm these predictions, C0i do 2 not. This result, together with the fact that the spreads of Fi are significantly lesser than those of Fi1 , suggests preference for the C0 estimates obtained by the method M2 here proposed. Two extreme γ values (0.80 and 0.40) , proposed by Angell et al. (1971) and Saffman (1963), were finally considered. The C0 values thus obtained lie in the range of values previously proposed in the literature by many authors. It is suggested that one cause of such different C0 values might be the different values of γ.

266

D. ANFOSSI ET AL.

Appendix A: Calculation of C0 According to M2 It is the aim of this appendix to derive in more details the M2 final formula (Equation (12) of the main text). Let us rewrite here Equations (8a), (10a), (9a) and (10b) of Section 2 (Monin and Yaglom, 1975), which are valid for the inertial subrange DiL (τ ) = C0i τ SiL (n) =

B0i −2 n 2π

(A1) n > nL

(A2)

DiE (τ ) = αi Cs ( u) ¯ 2/3 τ 2/3 SiE (n) =

(A3)

αi αu ( u) ¯ 2/3 n−5/3 (2π )2/3

n > nE ,

(A4)

where, i = u, v, w, DiL (τ ) and DiE (τ ) are the second-order Lagrangian and Eulerian velocity structure functions, SiL (n) and SiE (n) are the Lagrangian and Eulerian energy spectra, τ is the time lag, n is the frequency,  is the ensembleaverage rate of dissipation of turbulent kinetic energy, u¯ is the mean wind speed, αu = 1, αv = αw = 4/3 (Sorbjan, 1989) and B0i , C0i , Cs and αu are numerical constants. nL and nE may be regarded as the inverse of the characteristic time scale TiL and TiE . Following Corrsin (reference from Pasquill, 1974, p. 89) let us integrate Equations (A2) and (A4) from zero to infinite, obtaining Z ∞ B0i B0i  2 σLi = n−2 dn = , (A5)  2π nL 2π nL and 2 σEi

αi αu = ( u) ¯ 2/3 (2π )2/3

Z

∞

n nE

−5/3

3 αi αu dn = 2 (2π )2/3



 u¯ nE

2/3 .

(A6)

2 2 Since σLi and σEi are assumed equal (the turbulent kinetic energy is the same in both Lagrangian and Eulerian reference schemes), the elimination of  gives  3/2 nE (αi αu )3/2 u¯ TiL 3 = E = . (A7) nL 2 B0i σi Ti

Let us now derive expressions (8c) and (9c), originally proposed by Monin and Yaglom (1975). The starting point is given by the following general relationship relating, in the Lagrangian framework, the structure function, the autocorrelation function and the energy spectrum Z ∞ L L L Di (τ ) = 2[Ri (0) − Ri (τ )] = 2 [1 − cos(ωτ )]EiL (ω) dω, (A8) 0

267

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

where ω is the angular frequency. Changing from angular frequency ω to frequency n (by setting ω = 2π n and SiL (n) = 2π EiL (ω)), using Equation (A2) and integrating over n, gives Z B0i ∞ DiL (τ ) = [1 − cos(2π τ n)]n−2 dn = B0i π τ, (A9) π 0 where the following approximate result is used Z ∞ [1 − cos(x)] π dx ∼ = . 2 x 2 0

(A10)

Comparing Equations (A9) and (A1) we obtain B0i =

C0i , π

(A11)

that is exactly Equation (8c). With a similar procedure we can also obtain Equation (9c). By writing the general Eulerian relationship Z ∞ DiE (τ ) = 2[RiE (0) − RiE (τ )] = [1 − cos(ωτ )]EiE (ω) dω, (A12) 0

and changing from ω to n, one obtains Z ∞ 2αi αu 2/3 DiE (τ ) = ( u) ¯ [1 − cos(2π τ n)]n−5/3 dn = 4αi αu ( uτ ¯ )2/3 , (2π )2/3 0 (A13) where, in this case, the following approximate result is used Z ∞ [1 − cos(x)] dx ∼ = 2. x 5/3 0

(A14)

Comparing Equations (A13) and (A3) we obtain Cs αu ∼ , = 4

(A15)

that is exactly Equation (9c). By inserting Equations (A11) and (A15) into Equation (A7), we are now in a position to derive Equation (11) of Section 2, here named (A16), TiL = TiE

 3/2 3 (αi αu )3/2 u¯ π . 2 B0i σi

(A16)

268

D. ANFOSSI ET AL.

References Anand, M. S. and Pope, S. B.: 1985, ‘Diffusion behind a Line Source in a Grid Turbulence’, in L. J. S. Bradbury, F. Durst, B. E. Lauder, F. W. Schmidt, and J. H. Whitelaw (eds.), Turbulent Shear Flows, Vol. 4, Springer-Verlag, Berlin, pp. 46–61. Angell, J. K., Pack, D. H., Hoecker, W. H., and Delver, N.: 1971, ‘Lagrangian-Eulerian Time-Scale Estimated from Constant Volume Balloon Flights Past a Tall Tower’, Quart. J. Roy. Meteorol. Soc. 97, 87–92. Cassardo, C., Sacchetti, D., Morselli, M. G., Anfossi, D., Brusasca, G., and Longhetto, A.: 1995, ‘A Study of the Assessment of Air Temperature, and Sensible and Latent Heat Fluxes from Sonic Anemometer Observations’, Nuovo Cimento 18C, 419–440. Corrsin, S.: 1963, ‘Estimates of the Relations between Eulerian and Lagrangian Scales in Large Reynolds Number Turbulence’, J. Atmos. Sci. 20, 115–119. Csanady, G. T.: 1973, Turbulent Diffusion in the Environment, Geophysics and Astrophysics Monograps, Reidel, Boston, 248 pp. Degrazia, G. A. and Anfossi, D.: 1998, ‘Estimation of the Kolmogorov Constant C0 from Classical Statistical Diffusion Theory’, Atmos. Environ. 32, 3611–3614. Du, S., Sawford, B. L., Wilson, J. D., and Wilson, D. J.: 1995, ‘Estimation of the Kolmogorov Constant for the Lagrangian Structure Function, Using a Second-Order Lagrangian Model of Grid Turbulence’, Phys. Fluids 7, 3083–3090. Du, S., Wilson, J. D., and Yee, E.: 1994, ‘Probability Density Functions for Velocity in the Convective Boundary Layer and Implied Trajectory Models’, Atmos. Environ. 28, 1211–1217. Ferrero, E. and Anfossi, D.: 1998a, ‘Sensitivity Analysis of Lagrangian Stochastic Models for CBL with Different PDF’s and Turbulence Parameterizations’, in S. E. Gryning and N. Chaumerliac (eds.), Air Pollution Modelling and Its Applications XI, Vol. 22, Plenum Press, New York, pp. 267–273. Ferrero E. and Anfossi D.: 1998b, ‘Comparison of PDFs, Closures Schemes and Turbulence Parameterizations in Lagrangian Stochastic Models’, Int. J. Environ. Pollut. 9, 384–410. Flesch, T. K. and Wilson, D. J.: 1992, ‘A Two-Dimensional Trajectory Simulation Model for NonGaussian Inhomogeneous Turbulence within Plant Canopies’, Boundary-Layer Meteorol. 61, 349–374. Frisch, U.: 1995, Turbulence, Cambridge University Press, U.K., 296 pp. Garratt, J. R.: 1992, The Atmospheric Boundary Layer, Cambridge University Press, U.K., 316 pp. Gifford, F. A.: 1955, ‘A Simultaneous Lagrangian-Eulerian Turbulence Experiment’, Mon. Wea. Rev. 83, 293–301. Hanna, S. R.: 1981, ‘Lagrangian and Eulerian Time-Scale in the Daytime Boundary Layer’, J. Appl. Meteorol. 20, 242–249. Hay, J. S. and Pasquill, F.: 1959, ‘Diffusion from a Continuous Source in Relation to the Spectrum and Scale of Turbulence’, in F. N. Frenkiel and P. A. Sheppard (eds.), Atmospheric Diffusion and Air Pollution, Advances in Geophysics, Vol. 6, Academic Press, pp. 345–365. Hinze, J. O.: 1975, Turbulence, McGraw Hill, New York, 790 pp. Horst, T. W. and Weil, J.: 1992, ‘Footprint Estimation for Scalar Flux Measurements in the Atmospheric Surface Layer’, Boundary-Layer Meteorol. 59, 279–296. Hurley, P. J. and Physick, W.: 1991, ‘A Skewed Homogeneous Lagrangian Particle Model for Convective Conditions, Atmos. Environ. 25A, 1313–1325. Hurley, P. J. and Physick, W. L.: 1993, ‘A Lagrangian Particle Model of Fumigation by Breakdown of the Nocturnal Inversion’, Atmos. Environ. 27A, 619–642. Kaimal, J. C. and Finnigan, J. J.: 1994, ‘Atmospheric Boundary Layer Flows’, Oxford University Press, 289 pp. Kaiser, R. and Fedorovich, E.: 1998, ‘Turbulence Spectra and Dissipation Rates in a Wind Tunnel Model of the Atmospheric Convevtive Boundary Layer’, J. Atmos. Sci. 55, 580–594.

ESTIMATION OF THE LAGRANGIAN STRUCTURE FUNCTION CONSTANT C0

269

Landau, L. D. and Lifshitz, E. M.: 1987, Fluid Mechanics, Second Edition, Pergamon Press, Oxford, 539 pp. Luhar, A. K. and Britter, R. E.: 1989, ‘A Random Walk Model for Dispersion in Inhomogeneous Turbulence in a Convective Boundary Layer’, Atmos. Environ. 23, 1191–1924. McMillen, R. T.: 1988, ’An Eddy Correlation Technique with Extended Applicability to Non Simple Terrain’, Boundary-Layer Meteorol. 43, 231–245. Monin, A. S. and Yaglom, A. M.: 1975, Statistical Fluid Mechanics: Mechanics of Turbulence, Vol. 2, MIT Press, 874 pp. Moore, G. E., Liu, M. K., and Shi, L. H.: 1985, ‘Estimates of Integral Time Scales from a 100-m Meteorological Tower at a Plains Site’, Boundary-Layer Meteorol. 31, 349–368. Panofsky, H. A. and Dutton, J. A.:1984, Atmospheric Turbulence – Model and Methods for Engineering Applications, John Wiley and Sons, New York, 397 pp. Pasquill, F.: 1974, Atmospheric Diffusion, Wiley & Sons, New York, 429 pp. Physick, W. L., Noonan, J. A., McGregor, J. L., Hurley, P. J., Abbs, D. J., and Manins, P. C.: 1994, ‘LADM: A Lagrangian Atmospheric Dispersion Model’, Technical Report, 24, CSIRO Division of Atmospheric Research, Australia, Rodean, H. C.: 1991, ‘The Universal Constant for the Lagrangian Structure Function’, Phys. Fluids A3, 1479–1480. Rodean, H. C.: 1994, Notes on the Langevin Model for Turbulent Diffusion of “Marked” Particles, UCRL-ID-115869 Report of Lawrence Livermore National Laboratory, 122 pp. Rotach, M. W, Gryning, S. E., and Tassone, C.: 1996, ‘A Two-Dimensional Lagrangian Stochastic Dispersion Model for Daytime Conditions’, Quart. J. Roy. Meteorol. Soc. 122, 367–389. Saffman, P. G.: 1963, ‘An Approximate Calculation of the Lagrangian Autocorrelation Coefficient for Stationary Homogeneous Turbulence’, Appl. Sci. Res. 11, 245. Sawford, B. L.: 1991, ‘Reynolds Number Effects in Lagrangian Stochastic Models of Turbulent Dispersion’, Phys. Fluids A3, 1577–1566. Sawford, B. L.: 1993, ‘Recent Developments in the Lagrangian Stochastic Theory of Turbulent Dispersion’, Boundary-Layer Meteorol., 62, 197–215. Sawford, B. L. and Borgas, M. S.: 1994, ‘On the continuity of Stochastic Models for the Lagrangian Velocity in Turbulence’, Physica D 76, 297–311. Sawford, B. L. and Guest, F. M.: 1988, ‘Uniqueness and Universality of Lagrangian Stochastic Models of Turbulent Dispersion’, in 8th Symposium on Turbulence and Diffusion, Amer. Meteorol. Soc., San Diego, CA, pp. 96–99. Sawford, B. L. and Tivendale, C. M.: 1992, ‘Measurements of Concentrations Statistics Downstream of a Line Source in Grid Turbulence’, in Proceedings of the 11th Australasian Fluid Mechanics Conference, University of Tasmania, pp. 945–948. Slade, D. H.: 1968, Meteorology and Atomic Energy 1968, USAEC, Division of Technical Information Extension, Oak Ridge, U.S.A., 445 pp. Smith, F. B.: 1967, ‘Eulerian and Lagrangian Time-Scale Relationship in One-Dimensional Turbulence’, in USAEC Meteorological Information Meeting, Chalk River, Canada (AECL-2787), pp. 476–483. Sorbjan, Z.: 1989, Structure of the Atmospheric Boundary Layer, Prentice Hall, New Jersey, 317 pp. Tassone, C., Gryning, S. E., and Rotach, M.: 1994, ‘A Random Walk Model for Atmospheric Dispersion in the Daytime Boundary Layer’, in S. E. Gryning and M. Millan (eds.), Air Pollution Modeling and Its Application X, Plenum Press, pp. 243–251. Tennekes, H.: 1982, ‘Similarity Relations, Scaling Laws and Spectral Dynamics’, in F. T. M. Nieuwstadt and H. van Dop (eds.), Atmospheric Turbulence and Air Pollution Modelling, Reidel, Dordrecht, pp. 37–68. Thomsen, A., Barring, L., Gryning, S. E., Sogaard, H., and Thorgeirsson, H.: 1994, Data Report for the NOPEX Tisby Site – Concentrated Field Effort #1, May 27–June 23, 1994.

270

D. ANFOSSI ET AL.

Thomson, D. J.: 1987, ‘Criteria for the Selection of Stochastic Models of Particle Trajectories in Turbulent Flows’, J. Fluid Mech. 180, 529–556. Van Dop, H., Nieuwstadt, F. T. M., and Hunt, J. C. R.: 1985, ‘Random Walk Models for Particle Displacements in Inhomogeneous Unsteady Turbulent Flows’, Phys. Fluids 28, 1639–1653. Wandel, C. F. and Kofoed-Hansen, O.: 1962, ‘On the Eulerian-Lagrangian Transform in the Statistical Theory of Turbulence’, J. Geophys. Res. 76, 3089–3093. Weil, J. C.: 1990, ’A Diagnosis of the Asymmetry in Top-Down and Bottom-Up Diffusion Using a Lagrangian Stochastic Model’, J. Atmos. Sci. 47, 501–515. Wilson, J. D. and Flesch, T. K.: 1993, ‘Flow Boundaries in Random-Flight Dispersion Models: Enforcing the Well-Mixed Condition’, J. Appl. Meteorol. 32, 1695–1707. Wilson, J. D. and Sawford, B. L.: 1996, ‘Review of Lagrangian Stochastic Models for Trajectories in the Turbulent Atmosphere’, Boundary-Layer Meteorol. 78, 191–210. Wilson, J. D., Legg, B. J., and Thomson, D. J.: 1983, ‘Calculation of Particle Trajectories in the Presence of a Gradient in Turbulent-Velocity Variance’. Boundary-Layer Meteorol. 27, 163–169. Zannetti, P.: 1990, Air Pollution Modelling, Computational Mechanics Publications, 444 pp.