CHARASTERISTICS OF THE NEAR-SURFACE TURBULENCE DURING A BORA EVENT

Wind velocity at the town of Senj was measured at 13 m above the ground with a 3D ultrasonic anemometer operating at 4 Hz sampling frequency. The severe bora case that occurred on 07 January and lasted to 11 January 2006 is analyzed here. This data set is used for evaluation of the turbulent kinetic energy, TKE, and its dissipation rate, ε. Some considerations about defining turbulent perturbations of the bora wind speed are pointed out. The inertial dissipation method for estimation of ε is used. The empirical length scale parameter for this event is estimated with respect to ε and TKE.


Introduction
Downslope windstorms occur at different locations worldwide, as found in e.g.Ólafsson and Ágústsson (2007) in their brief review regarding this subject.They mention the Boulder windstorms in westerly flow in North America, windstorms in Norway in westerly flow, Greenland windstorms in both westerly and easterly flow, downslope windstorms in southerly flow in the Alps, as well as bora windstorms in northeasterly flow over the Dinaric Alps.Although these winds share some similar dynamical patterns, they also greatly differ due to the phenomena such as breaking gravity waves, hydraulic jumps, atmospheric rotors, gap jets, boundary-layer separation, flow splitting and non-stationary flow behaviour (e.g.Gohm et al., 2008).Due to possible severity these winds represent one of the biggest weather threats to vegetation, infrastructure, traffic and lives.Therefore, there is a great and continuous interest in the broader Correspondence to: Ž. Večenaj (zvecenaj@gfz.hr)meteorological community and significant amount of effort is invested in downslope windstorms studies.
Besides more pleasant airflows at the northeastern Adriatic coast, such as sea and land breezes (e.g.Prtenjak et al., 2008), there are two other, more vigorous wind types; sirocco (e.g.Pasarić et al., 2007) and bora (e.g.Jurčec, 1981).Bora (locally bura) is a downslope windstorm that blows from the northeastern quadrant in the lee of the coastal mountains when the relatively cold northeasterly flow impinges on the Dinaric Alps (e.g.Yoshino, 1976;Poje, 1992).The most frequent occurrence of bora is during the winter season with duration from several hours to several days (e.g.Jeromel et al., 2009).It possesses a wide spectrum of average wind speeds, and due to strong gustiness the speed maxima may surpass 60 m s −1 (e.g.Belušić and Klaić, 2006).
There is a large body of work on bora trying to explain its nature.The most recognized explanation can be found in Smith (1987).He pointed out the hydraulic behaviour of the mean bora flow.However, the turbulent characteristics of bora remain relatively unexplored.There are few papers in the literature regarding this subject.Belušić et al. (2006) analyzed bora from the data measured in the coastal town of Senj with a cup anemometer at 1 Hz sampling frequency.They found that the magnitude of high frequency variance (periods less than 1 min), representing the turbulence, is closely related to the mean wind speed magnitude.This indicates that the near-surface turbulence is primarily generated locally by friction and wind shear.Conversely, at lower frequencies, i.e. periods roughly between 1 and 10 min, a pulsating variability appears that is of non-local origin (Belušić et al., 2004(Belušić et al., , 2006(Belušić et al., , 2007)).
Turbulent kinetic energy (TKE) is one of the most important variables in micrometeorology because it is a measure of the intensity of turbulence.During a bora event, due to its dynamics and the consequent flow severity, the turbulence is strongly developed in the lee of the mountain.Belušić and Klaić (2006)  Individual terms in the TKE budget equation (e.g.Stull, 1988) describe physical processes that generate or suppress turbulence.One of those terms is the turbulence kinetic energy dissipation rate, ε, which describes dissipation of TKE by molecular viscosity into the heat.Information about ε, among the other, can provide an insight into the local turbulence strength.This finds a variety of applications, e.g.turbulence nowcasting at the airports (e.g.Frech, 2007).
Micro-scale properties of bora and other downslope windstorms, e.g.turbulence spectra and the consequent parameterization's refinement for fine-scale modelling, are still subjects of vigorous research (e.g.Baklanov and Grisogono, 2007;Grisogono and Belušić, 2008).Finding a proper set of mixing length-scales for parameterization of turbulent fluxes in numerical models under changes between weakly-and strongly-stratified boundary-layer conditions appears necessary (e.g.Mahrt, 1998;Mauritsen et al., 2007;Gohm et al., 2008;Zilitinkevich et al., 2008;Grisogono and Belušić, 2009).In particular, wave-turbulence interactions in and around breaking mountain waves (e.g.strong to severe bora), are sensitive to parameterizations and particular values of TKE dissipation (e.g.Epifanio and Qian, 2008).All these improvements needed for a more faithful turbulence parameterization in ever-refining models are impossible without a more complete knowledge of ε.This study deals with ε under strong bora conditions on the basis of the high frequency turbulence data which has been addressed in some details only recently (Belušić et al., 2006).A 4 day bora episode recorded in January 2006 is analysed here.A few preliminary results of this study can be found in Večenaj et al. (2007Večenaj et al. ( , 2009)).Until this dataset was produced, no such analysis could have been performed for the bora wind.
In the first part of the paper a theoretical background of the estimation of ε and its parameterization is pointed out.In the second part, ε is first estimated using inertial dissipation method (IDM) for the 4 h strongest bora interval within the 4 day episode.The significance of those estimations is evaluated by comparing ε and TKE over the entire bora episode.In addition, the dependence of ε on the mean streamwise velocity component, U , is investigated.Some considerations about the characteristics turbulent length scale, , are carried out.Finally, our concluding remarks are summarized in the last section.

Data
The 3-D wind velocity measurements were performed at Senj (44.99 • N, 14.90 • E, 2 m above the mean sea level) at a height of 13 m above ground with the WindMaster ultrasonic anemometer (Gill Instruments).The location of Senj is shown in Fig. 1.It is a small town at the northeastern Adriatic coast with buildings in general no higher than a few meters above the ground.The anemometer was mounted on the roof of the Senj sea-port captaincy building, several meters from the sea, one of the taller exceptions among the urban facilities.Concerning the urban influence, Roth (2000) states that the shapes of the velocity spectra over different urban sites are in good agreement with the reference rural data.Since our analysis is mostly based on power spectra, we find this to be a sufficient guarantee that our analysis is not significantly affected by the urban settings of Senj.
The anemometer recorded the data with a frequency of 4 Hz.Aliasing of the unresolved frequencies is a factor that might contaminate the spectra (e.g.Moore, 1986).Nevertheless, the true sampling frequency of the anemometer is 40 Hz and the data are block averaged to the outputting 4 Hz.This is so called "digital prefiltering", which successfully keeps aliasing under control (e.g.Kaimal and Finnigan, 1994).The observed bora episode extends from 7 to 11 January 2006 (a 4 day time series).During this period the ECMWF analysis shows the presence of a deep and persistent anticyclone above the NE Europe.This anticyclone was already formed on 7 January 2006 at 00:00 UTC and remained stationary up to 10 January 2006 at 00:00 UTC when the cyclone from the north started to repress it.Figure 2 depicts the 00:00 UTC MSL pressure and the 850, 700 and 500 hPa geopotential height fields on 9 January 2006.These synoptic conditions ensure a persistent supply of cold air impinging on the Dinaric Alps; thus the favourable conditions for the development of long-lasting downslope wind in the lee are provided (e.g.Jurčec, 1981;Heimann, 2001).
The mean wind direction during this event is ≈55 • .The coordinate system is rotated in the mean direction and the horizontal wind is decomposed into the along (streamwise) and cross (transverse) wind components.Figure 3a depicts the measured 4 day time series of the streamwise velocity component.The grey segment is the strongest 4 h interval in this event which is chosen for a detailed analysis (Fig. 3b). 1 0 0 5   1 0 1 0   1 0 1 5   10 15   1 0 1 5   1 0 1 5  1 0 2 0  1 0 2 0  1020   1 0 2  3 Theoretical background

Inertial dissipation method
In order to evaluate the TKE dissipation rate ε, IDM provided by the Kolmogorov's 1941 hypothesis can be employed (e.g.Tennekes and Lumely, 1972).It requires the validity of the Taylor's hypothesis on frozen turbulence (TH), which is needed for the transformation from time (frequency, f ) to space (wavenumber, k) domain.The criterion for the validity of TH is σ M < 0.5M, where M is the mean horizontal wind speed and σ M is the standard deviation (e.g.Stull, 1988).According to Kolmogorov, the power spectrum follows the −5/3 law in the inertial subrange: where the usual units are valid but neglected in Eq. (1) for simplicity (e.g.Albertson et al., 1997).The strongest 4 h interval (grey area in Fig. 3a) chosen for the detailed analysis with the 10 min mean superimposed.Using TH and Eq. ( 1), ε can be evaluated from (e.g.Champagne et al., 1977): where S u (f ) is the spectrum and α is the Kolmogorov constant for the velocity component.
According to Batchelor (1959) and Champagne (1978), a strong statement of local isotropy and the existence of the inertial subrange are given by the 4/3 ratios of spectra of the transverse (v) to streamwise (u) and the vertical (w) to streamwise (u) wind component.However, several authors (e.g.Champagne, 1978;Mestayer, 1982) show that the −5/3 law for the streamwise component spectrum can often be extended outside of the inertial subrange.This fact will be used as a justification for the data analysis presented here.

Parameterization of ε
In many atmospheric models "local closure" techniques are used for turbulence parameterization.Starting from the local one-and-a-half-order closure, ε can be parameterized using the mean value of TKE (e.g.Mellor and Yamada, 1974): where TKE= 1 2 u 2 + v 2 + w 2 and u , v and w are the turbulent perturbations of the streamwise, transverse and vertical velocity components, respectively, while overbars represent a suitable averaging.The parameter is the empirical length scale that is closely related to the size of dissipating turbulent eddies (a proportionality constant is absorbed in for convenience).By comparing TKE variations with those of ε and using Eq.(3), one can obtain an insight into the robustness and consistency of the ε estimation.

Analysis of the strongest bora interval
In order to evaluate ε, we used the IDM method discussed previously.This method assumes the existence of the inertial subrange which is characterized by local isotropy.Thus, a test of the local isotropy is performed on our bora wind data.The strongest bora 4 h interval of the 4 day time series is chosen (Fig. 2b) and the streamwise (S u ), the transverse (S v ) and the vertical (S w ) velocity spectra as functions of frequency are calculated using the fast Fourier transform (FFT).The spectra were smoothed by block averaging spectral amplitudes using 56 windows with 50% overlap.Each window contained 1024 data points.Figure 4a shows that the S v to S u ratio (grey marks) approximately reaches the theoretical value of 4/3 for f > 0.8 Hz, while S w to S u ratio (black marks) reaches 4/3 only at the end of the frequency band (which is the Nyquist frequency of 2 Hz).This suggests that an inertial subrange may exist within the frequency limits measured by the sonic anemometer 13 m above the ground.Piper and Lundquist (2004) found similar result in their data (shown in their Fig. 1). Figure 4b shows a loglog representation of S u (solid line), S v (dashed line) and S w (dashed-dotted line) from the 4 h interval.There is an indication of the −5/3 behaviour of spectra outside and nearby the inertial subrange, as revealed from the comparison with the −5/3 slope line.
Regardless of whether the inertial subrange exists within the current frequency limits (Fig. 4a), the IDM is applied to our data due to the extension of the −5/3 law discussed in Sect.3.1.The standard deviation of the total wind speed for the chosen 4 h interval is σ M =0.35 M, so the application of TH is justified.The Kolmogorov constant α=0.53 (e.g.Champagne, 1978;Oncley et al., 1996;Piper and Lundquist, 2004) is used in Eq. ( 2).Only streamwise velocity component perturbations are used for computing ε using IDM over the frequency band that best corresponds to the −5/3 slope.For this interval, U =15.1 m s −1 , ε=1.22 m 2 s −3 .

Analysis of the entire bora episode
The purpose of this section is to gain an insight into the evolution of TKE and ε in time as well as in their relationship.Also, we investigate the dependence of ε on U and estimate over the entire bora episode.In order to define the bora near-surface turbulent perturbations for the TKE calculation, it is important to find the scale that delimits the turbulence from the mean and/or mesoscale flow.This is not always straightforward, particularly in the complex flow such as the bora wind where several scales interact (e.g.Grubišić, 2004;Belušić et al., 2007;Gohm et al., 2008;Grisogono and Belušić, 2009).Hence, we first explore the proper value of the bora turbulence limit scale using different approaches.
The limit is usually defined in terms of spectra, more specifically, by an assumed energy gap at the mesoscale (e.g.Metzger and Holmes, 2008).However, the gap is not present in a number of situations; hence the difficulty in determining the scale.Figure 5 depicts the u component spectrum representing the mean value of spectra calculated from the 16 consecutive bora intervals within the 4 day time series that are 6 h long.Each spectrum was obtained from FFT and smoothed by block averaging spectral amplitudes using 3 windows with 50% overlap.A window contained 86 400 data points.While there is a noticeable energy gap at ∼10 min period, the energy peak appearing at slightly smaller periods (around 5 min) originates from the non-local quasi-periodic pulsations (Belušić et al., 2004).The pulsations are the manifestation of Kelvin-Helmholtz instabilities appearing in the approximately 1 km thick layer of strong wind shear above the bora jet, at heights above 500 m (Belušić et al., 2007).While these features are at the border between waves and turbulence, they are not related to local generating mechanisms and therefore should not be taken into consideration as the local turbulence.Since the peak related to pulsations is not sharp, it appears that the local turbulence clearly occurs only at periods ≤1 min, as in Belušić et al. (2006).The latter can be further substantiated as follows.
The relationship between the standard deviation of the streamwise velocity component, σ u , and U in the surface layer is usually linear (e.g.Stull, 1988), which was shown by Belušić et al. (2006) to be valid for the local bora turbulence.In order to examine the relationship for different limit scales, moving averages with lengths from 1 min to 1 h have been subtracted from the entire 4 day bora episode to determine the perturbations.Then a power law of the form: was fitted to the scatter diagrams of σ u vs. U for different moving averages (the fitting procedure is described near the end of this section).Figure 6a depicts the scatter diagram with σ u obtained by subtracting the 1 min moving average, and the fit whose power coefficient b is 1.025.The values of σ u and U were averaged for each hour of the whole 4 day bora episode.Figure 6b shows the power coefficient b relative to the moving average length.Expecting b to be 1, it is clear that the closest value is achieved for the 1 min moving average.Therefore, based on the two approaches above we have decided to use the 1 min moving average for determination of the local turbulence perturbations in this bora episode.Due to similarity with other bora episodes studied by Belušić et al. (2004Belušić et al. ( , 2006)), it seems that this limit scale may be recommended for other bora cases.
The IDM method for the ε evaluation presented in Sect. 3 is applied to the entire bora episode.The 1 min running average is subtracted from the observed bora episode and the residuals are defined as turbulent perturbations.This newly formed bora time series is divided into 96 intervals of 1 h length.For each of the intervals, the possibility of TH implementation is tested by calculating the criterion required by TH.All the intervals satisfy this criterion.As before, spectra are calculated using FFT, but this time smoothed by block averaging of 14 windows with 50% overlap.Each window contains 1024 data points.
Figure 7a depicts the time evolution of TKE, and ε evaluated independently.It is seen that the time series closely follow each other, which indicates a significant degree of robustness in the ε estimation.The TKE values in the strongest part of the bora episode reach 20 J kg −1 .Figure 7b  time evolution of estimated from Eq. ( 3).The mean value of during the entire bora episode is ≈60 m, which indicates the average eddy size at which ε is most effective in dissipating TKE at the height of the sensor.
Since TKE is defined as in Eq. ( 3), and σ u ∝ U , one can easily show that TKE∝ U 2 .Therefore, from Eq. ( 3) it follows that ε ∝ U 3 .All these relations can be written as: where a and b are coefficients while x and y represent a certain pair of variables.Here we also define theoretical values of b as b t =2, for x = U and y=TKE; and b t =3, for x = U and y = ε.Likewise, for Eq. ( 3) by definition b t =1.5.These relations (i.e., coefficients) are tested by fitting scatter diagrams of ε vs. TKE, TKE vs. U , and ε vs. U , using the least squares method.In order to evaluate properties of an estimator (Eq.4b), 100 different samples containing 80 % randomly chosen data points from a certain scatter diagram are constructed (the bootstrapping method).Each of these samples is fitted to obtain the values of coefficients a and b for the estimator.Finally, we apply the linear regression to calculate the correlation coefficient squared, r 2 , which provides the efficacy of the fit.For a comparison with the previous theoretical values of the coefficients b t , the whole procedure is repeated, now with a priori assigned power coefficients b t of the estimator: The corresponding correlation coefficient squared for Eq. ( 4c) is r 2 t .The results are shown in Table 1 and Fig. 8. Figure 8 depicts the scatter diagrams.The solid line denotes the fit with a priori unknown coefficients, while the dashed line denotes the fit with the theoretically assigned coefficient b t .The correlation coefficient for the ε vs. TKE scatter diagram (Fig. 8a) is relatively high for both fits (r 2 =0.912 for the a priori unknown coefficients a and b; r 2 t =0.903 for the theoretically assigned coefficient b t ), which indicates a strong relation between the two variables.While theory requires for ε to be proportional to TKE 1.5  as in Eq. (3), our data suggest the ε being proportional to ≈TKE 1.3 .The length-scale parameter , calculated as (a t ) −1 , where a t is the linear coefficient of the dashed fit, is 60.2 m.The correlation coefficients for TKE vs. U (Fig. 8b) and ε vs. U (Fig. 8c) for both fits are not as high, but are still significant.
In order to compare our results for ε with the similarity estimate, we have calculated u * and used the values to normalize ε as defined e.g. in Piper and Lundquist (2004): where u 2 * = u w 2 + v w 2 , k is the von Karman's constant, taken here to be 0.4 (e.g.Albertson et al., 1997), and h is the height of the instrument above ground.Since we did not measure temperature, we can not estimate stability range for this bora episode.However, it is reasonable to assume neutral to weakly stable surface layer during a strong bora event far down the mountain lee (e.g.Grisogono and Belušić, 2009).Wyngaard and Coté (1971) expect for ε in the neutral atmospheric conditions to be 1.We have fitted the scatter diagram of ε vs. u * (Table 1; Fig. 8d) using Eq. ( 4c) with the a priori assigned coefficient b t =3 and calculated the normalized ε as  Wyngaard and Coté (1971) and in agreement with the values from Piper and Lundquist (2004), their Fig. 8.

Comparison with parameterizations in numerical models
Many current state-of-the-art atmospheric numerical models still use the Blackadar length-scale parameterization, defined as in e.g.Mellor and Yamada (1974): where the constant C ≈ 5.3 for the neutral boundary layer, qdz, the constant β=0.1, z is the height above ground and q = √ 2TKE.We have compiled a number of vertical TKE profiles from different bora episodes and from two different models: MEMO6 (Mesoscale Model, version 6) mesoscale model, previously used for simulations of several bora events (Belušić and Klaić, 2004), and the WRF-ARW (Weather Research and Forecasting -Advanced Research WRF) model using the Mellor-Yamada-Janjic turbulence parameterization scheme (Večenaj et al., 2009).
The values of were calculated on the basis of these profiles.For each of the considered cases it amounts to ≈25 m.The model-based obviously underestimates the one derived here from TKE and ε values.This may imply an inadequacy of the Blackadar type of parameterizations for the bora related turbulence.

Discussion and conclusions
Certain properties of the near-surface turbulence related to the bora downslope windstorm are addressed.For the first time, the TKE dissipation rate, ε, during a bora event is estimated reliably.A typical bora case lasting for a couple of days in the town of Senj that is well-known for the frequent bora occurrence is presented.The analysis is performed using 4 Hz wind data.
Following Kaimal and Finnigan (1994), the ideal choice of sampling frequency is the one that reveals at least two octaves of the inertial subrange.With respect to h and U in our measurement setup, that would be ≈10 Hz.Thus, for a more detailed study of the bora inertial subrange and more conclusive results about the bora related ε, measurements of bora wind with at least 10 Hz sampling frequency are needed in future.
Although there was an initial doubt about the existence of the isotropic inertial subrange within frequencies available from the data, it is seen that the −5/3 slope extends sufficiently towards lower frequencies.Estimations of ε using the inertial dissipation method, IDM, applied to the streamwise velocity component and normalized with u 3 * , agree well with previous studies.Application of IDM to the transverse or vertical velocity component using a modified Kolmogorov constant, α m =4/3 α (e.g.Oncley et al., 1996), shows similar values of ε to those derived from the streamwise component.
The TKE values during the bora episode addressed here reach 20 J kg −1 .This is comparable with the values from Belušić and Klaić (2006).The TKE time series closely follow the IDM values of ε, which indicates the robustness of the ε estimate.
As expected, both ε and TKE increase with increasing the mean streamwise velocity component, U .Our data suggest that in this bora episode ε is proportional to ≈TKE 1.3 .This model explains 91% of ε vs. TKE variance.The Eq. (4c) fit to the ε vs. TKE scatter diagram with the a priori proposed power coefficient b t =1.5, as suggested in Eq. ( 3), leads to the length-scale parameter of ≈60 m.This value appears to be larger than those obtained from atmospheric numerical models, which indicates a possibility for improvements of turbulence parameterizations in models.
Future work related to the near-surface bora turbulence will include the analysis of the 4 Hz data for a variety of bora episodes, categorized by their nature (type), strength and seasonal period.In this way, a comprehensive picture of the bora-related near-surface turbulence will be revealed.At the same time, new field campaigns with sufficiently high sampling rates should be planned for a more thorough insight into the nature of bora turbulence.

Fig. 1 .
Fig. 1.The position of the town of Senj (indicated by black square) where the 4 Hz bora wind measurements were performed.The terrain height is given every 300 m.

Fig. 3 .
Fig. 3. (a) A 4 day raw 4 Hz data time series (from 7 to 11 January 2006) of the streamwise wind component measured in Senj, with the 1 h mean superimposed.(b)The strongest 4 h interval (grey area in Fig.3a) chosen for the detailed analysis with the 10 min mean superimposed.
Fig. 4. (a)The ratio between the transverse v and streamwise u velocity spectra (grey marks) and between the vertical w and streamwise u velocity spectra (black marks) showing the approach to the 4/3 ratio required by local isotropy (horizontal solid line), and (b) a log-log representation of u (thin solid line), v (thin dashed line) and w (thin dashed-dotted line) velocity power spectrum densities for the strongest 4 h bora interval.The black thick dashed line is the −5/3 slope.The v and w spectra have been reduced by a factor of 4 and 10, respectively, for presentation.

Fig. 5 .
Fig. 5.A lin-log representation of the mean power spectrum of the streamwise velocity component (solid curve) for the entire bora episode.The vertical lines denote the periods of 10 min (dotted line), 5 min (dashed-dotted line) and 1 min (dashed line).Vertical error bar shows the 95% confidence interval.
Fig. 6.(a) The scatter plot of σ u , obtained by subtracting the 1 min moving average from the entire bora episode, vs. U .The solid line represents the Eq.(4a) fit.(b) The power coefficients b from (a), for different moving averages subtracted from the entire bora episode.
Fig. 7. (a) Time evolution of TKE (black marks) and ε evaluated by IDM (grey marks), and (b) time evolution of the empirical length scale parameter estimated from Eq. (3) during the entire 4 day bora episode.

Fig. 8 .
Fig. 8. Scatter plots of (a) ε vs. TKE, (b) TKE vs. U , (c) ε vs. U and (d) ε vs. u * for the entire bora episode.Solid lines are fits with the a priori unknown coefficients a and b from the relation (4b).Dashed lines are fits with the a priori known coefficient b t from the relation (4c).

Table 1 .
Results of the fitting procedure.Each row contains fitting coefficients for a certain pair of variables addressed in the first column.These variables are evaluated over the 4 day bora episode.Given are the mean values and standard deviations of the a priori unknown coefficients a and b from the relation (4b), and the coefficient a t from the relation (4c), where b t is the a priori known theoretical coefficient.Also, the corresponding correlation coefficients squared r 2 and r 2 t are shown.= a t kh = 0.93.This is reasonably close to the expected value of ε