atmosphere - MDPI

atmosphere Article

Micro-Scale Properties of Different Bora Types Vinko Šoljan 1, * ID , Andreina Belušić 2 ID , Kristina Šarović 3 , Irena Nimac 4 , Stjepana Brzaj 4 , Jurica Suhin 1 , Martin Belavić 4 , Željko Veˇcenaj 2 and Branko Grisogono 2 1 2 3 4

*

Croatia Control Ltd., 10410 Velika Gorica, Croatia; [email protected] Department of Geophysics, Faculty of Science, University of Zagreb, 10000 Zagreb, Croatia; [email protected] (A.B.); [email protected] (Ž.V.); [email protected] (B.G.) EKONERG—Energy and Environmental Protection Institute, 10000 Zagreb, Croatia; [email protected] Meteorological and Hydrological Service, 10000 Zagreb, Croatia; [email protected] (I.N.); [email protected] (S.B.); [email protected] (M.B.) Correspondence: [email protected]; Tel.: +385-98-891-555

Received: 24 January 2018; Accepted: 16 March 2018; Published: 21 March 2018

Abstract: In this paper we use 20 Hz wind measurements on three levels (2, 5, and 10 m) to investigate the differences in micro-scale properties of different bora types, i.e., deep and shallow bora with further subdivision to cyclonic and anticyclonic bora cases. Using Fourier spectral analysis, we investigate a suitable turbulence averaging scale and bora gust pulsations. The obtained data set is further used to test the Monin–Obukhov similarity theory, the surface layer stratification, the behavior of the terms in the prognostic turbulence kinetic energy equation, and the wind profiles. One of our main goals is to identify possible micro-scale differences between shallow and deep bora types because of the possible different mountain wave dynamics in those flows. We found that a turbulence averaging scale of 30 min is suitable for this location and is in agreement with previous bora studies. The wind speed power spectral densities of all selected bora episodes showed pulsations with periods of 2–8 min. This suggests that mountain wave breaking was present in all cases, regardless of flow depth and synoptic type. The stability parameter analysis confirmed the near-neutral thermal stratification of bora; a consequence of intensive mechanical mixing. No significant differences related to bora type were observed in other micro-scale parameters. Keywords: bora; turbulence; micro-scale; pulsations; wind profiles; surface layer; mountain wave

1. Introduction Bora is a gusty downslope wind blowing from the northeast in the lee of the Dinaric Alps and other dynamically similar parts of the world [1]. It is dynamically generated by the interaction of airflow and orography [2]. Bora macro-scale characteristics have been investigated since the middle of the 20th century [3]. Early mesoscale research focused on the katabatic model of bora [4], which was later shown to be deficient in explaining stronger bora events. A major breakthrough in bora mesoscale research started with the ALPEX (Alpine Experiment) project in 1981 [5] and the subsequent findings by Smith [2], which showed that bora is essentially a dynamically generated wind, best explained by the hydraulic and wave breaking theory [6,7]. In recent years, bora’s mesoscale characteristics have been extensively studied during the Mesoscale Alpine Programme (MAP) project [8–10]. A recent review of bora mesoscale properties at the northeastern Adriatic can be found in [11]. With respect to the synoptic setup, three types of bora have been identified in the past: cyclonic, anticyclonic [3], and frontal [12]. The typical setup for cyclonic bora is when a mid-latitude cyclone moves to the southern Adriatic, pulling colder air from the continent to the eastern Adriatic coast. Anticyclonic bora blows when there is an anticyclone situated over Central Europe, extending over the Atmosphere 2018, 9, 116; doi:10.3390/atmos9040116

www.mdpi.com/journal/atmosphere

Atmosphere 2018, 9, 116

2 of 25

Dinaric Alps [3]. Frontal bora is characterized by a sudden increase in wind speed and short duration, following the passage of a cold front [12]. Both cyclonic and anticyclonic bora can be either deep or shallow, depending on the depth of synoptic flow over the mountains (e.g., [13]). For example, a common feature above a mature mid-latitude cyclone is the southwesterly flow in the divergent (eastern) flank of the upper-level trough (i.e., the flow from positive vorticity maximum to negative vorticity maximum at 300 hPa). This is why most cyclonic boras are shallow [2,4,14,15]. However, in the case of an occluded cyclone, there is usually a cut-off low in the upper levels, sometimes favorably aligned with the surface cyclone, thus providing deep NE (northeasterly) to N (northerly) flow throughout the troposphere. Anticyclonic deep bora occurs when a deep positively tilted trough passes and the upper N or NE flow in its western flank is above the southeastern quadrant of a surface anticyclone. Alternatively, the northwestern quadrant of a cut-off low can also provide deep NE flow above the surface anticyclonic bora. Anticyclonic deep bora seems to be the most frequent among deep bora types [16] because its synoptic setup is more common than that of cyclonic deep bora. Since strong anticyclones can persist for days, it is not uncommon for this type of bora to last up to a week. Deep bora is associated with vertically propagating mountain waves [2], while shallow bora is associated with wave breaking and violent downslope windstorms [7]. Shallow bora does not allow significant vertical propagation of wave energy, thus generating strong downslope windstorms in the lee. Vertical wind shear plays a significant role in the vertical propagation of waves in deep bora. In the case of positive wind shear (wind speed increasing with height), wave breaking does not occur at least until tropopause, because of the linearizing effect of the increasing wind speed [11]. In the case of weak vertical shear or wind speed decrease with height, wave breaking is likely to happen in the lower or middle troposphere, again reflecting mountain wave energy to lower levels and generating violent downslope windstorm. Regardless of more than a century of intensive research of bora climatology (e.g., [17–19]) and bora macro- and mesoscale properties (e.g., [11,20–22]), some important details about bora micro-scale properties are not yet known. One of them includes detailed characteristics of severe bora episodes. The most severe bora episodes (downslope windstorms), with gusts reaching up to 70 m·s−1 , are caused by wave breaking when there is a critical level above the mountaintop [7]. The critical level is usually marked by strong inversion [23] and a decrease in wind speed or change in wind direction by height, thus acting as an efficient reflector of wave energy. The critical level can be imposed by synoptic scale or generated by wave breaking itself—caused by wave amplitude increasing with height [24]. Severe bora typically induces shooting flow in the lee of coastal mountains [24] that may extend out over the sea in the form of multiple low-level jets behind mountain passes, while lee wakes (weaker flow regions) occur behind mountaintops [9]. Sea surface SAR (Synthetic Aperture Radar) data analysis by Kuzmić et al. [13] revealed the existence of secondary bora jets—caused by smaller mountain and island features (gaps and flanks)—that are only a few kilometers apart and several kilometers long. Moreover, they documented fine-scale convective cells pertaining to cold bora outbreak over relatively warm sea. Bora pulsations with periods of 3–4 min were first mentioned in the work of Watanabe [25], based on the experience of local fishermen. The first confirmation of those observations in the measured data was in the work of Petkovšek [26,27], who found pulsations with periods between 3 and 11 min. Although the existence of pulsations has been known for a long time, the detailed physics behind the gustiness and pulsations of bora has been addressed only recently. Belušić et al. [28] also found that the pulsations occur with periods between 3 and 11 min in the town of Senj, a location well known as a bora maximum site (Figure 1). Furthermore, this was also confirmed by using fine-scale numerical modeling [29], and measurements at the Pometeno Brdo (in a free translation, Pometeno Brdo means Swept Away Hill) [30]—a bora site upwind of the city of Split (Figure 1) that is about 200 km southeast from Senj. The former study found that the generation of gust pulsations was associated with mountain wave breaking and Kelvin-Helmholtz shear instability (KHI) above

Atmosphere 2018, 9, x FOR PEER REVIEW Atmosphere 2018, 9, 116

3 of 25 3 of 25

gust pulsations was associated with mountain wave breaking and Kelvin-Helmholtz shear instability (KHI) above bora shooting was flow.first Thisdemonstrated mechanism was first demonstrated by[31] Peltier the bora shooting flow. the This mechanism by Peltier and Scinocca for and Scinocca [31] for in mountain windstorms in Boulder, and CO,numerical USA. Measurements and numerical mountain windstorms Boulder, CO, USA. Measurements modeling studies [28,29,32] modeling studies [28,29,32] alsodisappear showed that thepresence pulsations disappear in thewind presence positive also showed that the pulsations in the of positive vertical shearof above the vertical wind shear above the mountaintop (e.g., the presence of an upper-tropospheric jet stream). mountaintop (e.g., the presence of an upper-tropospheric jet stream).

(a)

(b)

(c) Figure Figure 1.1. (a) (a)The TheAdriatic Adriatic coast coast with withlocations locations of of previous previous bora bora studies studies marked marked (Senj (Senj and and Pometeno Pometeno Brdo city of ofSplit). Split).The The Zadar and Zagreb Maksimir sounding stations are also marked. Brdo near the city Zadar and Zagreb Maksimir sounding stations are also marked. The The represents area of this study. shows thezoomed zoomedininview viewof ofthat thatarea, area, while while (c) shows boxbox represents thethe area of this study. (b)(b) shows the shows the the zoomed in view of the measurement site at the Maslenica Bridge (shown as the box in b). The Velebit zoomed in view of the measurement site at the Maslenica Bridge (shown as the box in b). Velebit Mountain Mountain is is mostly mostly north northof ofthe thesite. site.

Micro-scale Micro-scale characteristics characteristics of of severe severe turbulence turbulence in in the the wave wave breaking breaking region region are are the the focal focal point point of current bora research. In order to improve turbulence parameterization schemes in numerical of current bora research. In order to improve turbulence parameterization schemes in numerical models, cenaj et models, Veˇ Večenaj et al. al. [33] [33] evaluated evaluated turbulence turbulence kinetic kinetic energy energy (TKE) (TKE) and and its its dissipation dissipation rate rate for for aa bora cenaj et bora event event in in the the town town of of Senj. Senj. Veˇ Večenaj et al. al. [34] [34] estimated estimated the the turbulence turbulence dissipation dissipation rate rate along along the the Adriatic sea coast, using 4 Hz aircraft and dropsonde data obtained during the MAP project. For the Adriatic sea coast, using 4 Hz aircraft and dropsonde data obtained during the MAP project. For the Pometeno [35] analyzed bora wind speed profiles from 5 Hz5 data and found that PometenoBrdo Brdosite, site,Lepri Leprietetal.al. [35] analyzed bora wind speed profiles from Hz data and found they agreed with with commonly used used empirical power-law and the profiles. TheyThey also that they agreed commonly empirical power-law andlogarithmic-law the logarithmic-law profiles. found that thermal stratification of the surface layer is near neutral due to strong mechanical mixing. also found that thermal stratification of the surface layer is near neutral due to strong mechanical Using theUsing same the data, Lepri et al.Lepri [1,36]etfurther investigated turbulence intensity, Reynolds shear stress mixing. same data, al. [1,36] further investigated turbulence intensity, Reynolds and turbulence scale profiles for the mentioned shear stress andlength turbulence length scale profiles for thelocation. mentioned location. Without such high-frequency in situ measurements of wind in (e.g., spaceaircraft (e.g., Without such high-frequency in situ measurements of wind speedspeed in space aircraft measurements) in time (single-point ground-based measurementson, on,e.g., e.g., meteorological meteorological measurements) and inand time (single-point ground-based measurements towers/masts), the exploration of bora micro-scale properties would not be possible. towers/masts), the exploration of bora micro-scale properties would not be possible. For For aa more more comprehensive insight into the nature of bora turbulence, even higher sampling frequency comprehensive insight into the nature of bora turbulence, even higher sampling frequency (e.g.,(e.g., >10 >10 measurements needed. hints at the of this study. Hz)Hz) measurements are are needed. ThisThis alsoalso hints at the goalgoal of this study.

Bora has a major influence on all forms of transportation, engineering structures, electrical and telecommunication grids, agriculture, sea dynamics, air pollution, tourism, and firefighting.


4 of 25

Bora has a major influence on all forms of transportation, engineering structures, electrical and telecommunication grids, agriculture, sea dynamics, air pollution, tourism, and firefighting. Engineering structures in areas prone to severe downslope windstorms must be strong enough to withstand these hurricane force winds. Agriculture in those areas must also be adapted to such harsh conditions. Transportation is the most vulnerable human activity, since severe bora episodes can completely shut down all road traffic to and from the coast. In some extreme cases, even the air traffic at the whole eastern Adriatic coast can be completely suspended. The Maslenica Bridge is a very important transportation route, connecting the southern and central Croatian coast—the northeastern Adriatic coast—with inland parts of Croatia. The purpose of this study is to test whether some of the previous results obtained for different measuring sites apply to the Maslenica Bridge location. Namely, the turbulence averaging time scale, bora pulsations, thermal stratification, TKE budget, and wind speed profiles. Furthermore, we aim to identify possible differences in those micro-scale properties of different bora types. As this has not been attempted before, it could give new insights into the turbulence characteristics of bora wind. For this purpose, we classify bora episodes by the flow depth and the synoptic type. As already mentioned, the flow depth is important in defining the mountain wave dynamics. We think that the synoptic type (i.e., cyclonic, anticyclonic or frontal) can influence the micro-scale properties of bora mainly because of different wind speeds, but also with different vertical wind and temperature profiles (e.g., stronger inversions in anticyclones), which additionally define mountain wave dynamics. The maximum wind speeds depend on the synoptic type because anticyclones have horizontal pressure gradient limit (inertial instability) and thus can never have wind speeds as high as very deep cyclones. Finally, the flow depth itself is also dependent on synoptic situation. In the following sections, we will explain all the methods and data used, show and discuss the obtained results, summarize the main findings, and provide conclusions. 2. Methods and Data 2.1. Measurement Site and Instruments The measurement site (15.53◦ E, 44.24◦ N, 78 m ASL) is settled ≈30 km northeast of the city of Zadar and ≈200 m northeast from the Maslenica Bridge on the A1 section of the Croatian motorway (Figure 1). High-frequency data were collected on a 10 m mast. WindMaster Pro ultrasonic anemometers (Gill Instruments Ltd., Lymington, UK) measured the 3D wind speed and sonic temperature at 2, 5, and 10 m levels above the ground during the period from 8 October 2015 to 11 February 2016. The data were sampled with a frequency of 20 Hz. This is the highest sampling rate for any prolonged bora in situ measurements that exist today. The anemometers were connected to a CR3000 data logger (Campbell Scientific Inc., Logan, UT, USA) and the whole system was powered by two 60 W solar panels. The ground in the immediate vicinity of the mast is characteristic for the Adriatic coastline, with prevailing bare rocks and some herb cover in the form of garrigue (low, shriveled, light scrub) and maquis (dense hard-leaf shirk). 2.2. Data Quality Check The measured data were quality checked and despiked using two methods described in [37]. Unrealistic data values were detected by using absolute limits (100 m·s−1 for the wind speed and 100 ◦ C for the temperature) and linearly interpolated. Spikes were identified as three or fewer consecutive points in the time series with an amplitude larger than 3.5 standard deviations from the moving mean with a window length of 6000 data points (5 min). Spikes were also linearly interpolated and the process was repeated after increasing the standard deviation factor by 0.1 until no more spikes were detected. While interpolating the removed unrealistic data values and spikes, we kept a record of where those missing data points were located. We did this because we noticed that the interpolated continuous blocks of missing data had a negative influence on the spectral analysis (random missing data did not


5 of 25

have such large influence). After despiking, we downsampled the data to 10 Hz by averaging two Atmosphere 2018, 9, x FOR PEER REVIEW 5 of 25 consecutive data points in order to reduce the number of missing data points. After these procedures, the missing data were episodes, reducedwith to less than 1% of forone allepisode analyzed episodes, with of one 1% for all analyzed an exception where there was 3.9%an ofexception missing data episode where there was 3.9% of missing data at the 5 m level. at the 5 m level. 2.3. Criteria for Bora Episode Detection 2.3. Criteria for Bora Episode Detectionand andSelection Selection To identify borabora episodes ininthe we used used10-min 10-min averages u and v wind To identify episodes therecorded recorded data, data, we averages of uofand v wind components the m level. The direction wind direction distributions of different wind speed components fromfrom the 10 m 10 level. The wind distributions of different wind speed categories categories visualized as the wind2)rose (Figure 2) clearly thewind dominant wind directions of bora visualized as the wind rose (Figure clearly indicate the indicate dominant directions of stronger stronger bora events during the measurements. events during the measurements.

Figure 2. The wind roserose from 10-min of all alldata. data.The The wind speed category limits Figure 2. The wind from 10-minwind windaverages averages of wind speed category limits are inare in 1 , and m·s−m· thethe numbers ononthe therelative relativefrequency frequency of occurrence. s−1, and numbers theplot plotindicate indicate the of occurrence. −1from ◦ –60◦ .weThus, highest relative frequencyof of wind wind speed s−1·sis directions 15°–60°. The The highest relative frequency speed>10 >10m·m is from directions 15Thus, ◦ ◦ decided to to set set the criterion to 40° ± 40°. wind speed waslimit set towas 4 set ± The we decided the wind winddirection direction criterion to 40 40 .lowest The lowest windlimit speed −1 1in order to filter out weaker katabatic or other thermally driven flows, while still being able to m· to 4 m ·ss− in order to filter out weaker katabatic or other thermally driven flows, while still being −1) was used by Lepri et al. [35]. catch certain feeble onsets of bora. A similar threshold (5 m· able to catch certain feeble onsets of bora. A similar threshold (5s m ·s−1 ) was used by Lepri et al. [35]. Furthermore, the detection algorithm required that these conditions must be satisfied for at least 6 Furthermore, the detection algorithm required that these conditions must be satisfied for at least h, with an allowed discontinuity of 1 h. This is to allow for possible weaker periods or even wind 6 h, with an allowed discontinuity of 1 h. This is to allow for possible weaker periods or even wind reversals within the bora episodes caused by lee rotor formation or low-level flow separation reversals withinBora the episodes bora episodes caused leewere rotoralso formation low-level separation [22,38,39]. detected in thisby way visuallyor checked (not flow shown). We also [22,38,39]. tried Bora varying episodesthe detected in this way were also visually checked (not shown). We also varying detection settings, which confirmed that the stated settings were optimal tried because they the detection settings, which confirmed the stated were optimal because they caused minimal caused minimal fragmentation of that seemingly wholesettings bora episodes. Usingof these criteria, whole a total of 14 bora episodes were detected. For each of these episodes, the fragmentation seemingly bora episodes. synopticthese situation was aanalyzed surface analysiswere [40], detected. 500 hPa geopotential andthese meanepisodes, sea Using criteria, total ofusing 14 bora episodes For each of level pressure analysis [41] (NCEP GDAS/FNL—National Centers for Environmental Prediction the synoptic situation was analyzed using surface analysis [40], 500 hPa geopotential and mean Global Data Assimilation 0.25° × 0.25° analysis), and soundings from ZadarPrediction and sea level pressure analysis [41]System/Final (NCEP GDAS/FNL—National Centers for Environmental Zagreb stations [42] (University of Wyoming database). Figure 1 shows the position of the sounding Global Data Assimilation System/Final 0.25◦ × 0.25◦ analysis), and soundings from Zadar and Zagreb stations in relation to the measurement site. According to this analysis, we classified all bora stations [42] (University of Wyoming database). Figure 1 shows the position of the sounding stations episodes and selected the ones that unambiguously fell into one of the main bora type categories in relation (Table to 1). the measurement site. According to this analysis, we classified all bora episodes and selected the ones that unambiguously fell into one of the main bora type categories (Table 1).

Table 1. The main bora type categories.

Table 1. The main bora type categories. Bora Type Abbreviation shallow cyclonic SC Bora Type Abbreviation shallow anticyclonic SA shallow deep cyclonic cyclonic DCSC shallow anticyclonic deep anticyclonic DASA deep cyclonic DC frontal F deep anticyclonic DA frontal F


6 of 25

To emphasize the relevance of the flow depth in defining the mountain wave and bora dynamics, in this work, we propose a classification of bora episodes first by flow depth and then by its synoptic setup. The criterion for determining the depth of bora flow from sounding data was that the upper wind direction must be from the direction 40◦ ± 45◦ , which is similar to the one used by Smith [2]. An episode was classified as deep if this criterion was satisfied at least up to 500 hPa. Since the sounding data were available only at 00 UTC (Coordinated Universal Time) and 12 UTC, a 500 hPa geopotential height analysis was also used to subjectively assess the flow depth. Many episodes are transitional, and so may include change from cyclonic to anticyclonic (e.g., SC to SA; see Table 1). In some cases, this is also accompanied by change of the flow depth (e.g., SC to DA). In transitional episodes with one dominant type, that type is emphasized with bold font (see Table 2 in the Results and Discussion Section). 2.4. Bora Turbulence Spectra Before analyzing the turbulent characteristics of the selected bora episodes, an appropriate averaging time scale needs to be determined. The averaging time scale is needed to separate turbulent perturbations from the mesoscale and synoptic atmospheric motions. To separate the signals into mean and fluctuating components, we used Reynolds decomposition. The basic assumption of Reynolds decomposition is the existence of a local minimum (spectral “gap”) in the power spectral densities of various turbulence quantities, such as wind speed or potential temperature [43,44]. Before Reynolds decomposition, the mean wind direction was determined for each bora episode and the coordinate system was rotated, so the x-axis points downstream of this mean wind speed. The power spectral densities were calculated for the horizontal components (u and v) of wind speed using the Welch algorithm [45] and then smoothed using the frequency window that expands in width with frequency [46]. This frequency smoothing is needed to obtain representative spectral curve from the estimates, which exhibit excessive crowding and large scatter at the high-frequency end on a logarithmic scale. Since datasets of bora episodes contain blocks of missing data, the calculation of power spectral densities was not simple. For blocks of missing data shorter than 1 s (10 data points), the missing data were replaced by linearly interpolated values. If an episode contained blocks of missing data larger than 1 s, the episode was split into smaller segments with a minimum length of 3 h and without missing data. Then the power spectral densities were calculated for every segment using a window length equal to half the length of the smallest segment of one bora episode. Finally, a spectrum for a single bora episode was created by averaging spectral estimates of every segment in the episode. 2.5. Taylor Hypothesis After finding a suitable averaging period by spectral analysis, the time series was divided into non-overlapping block intervals with the length of the mentioned period. All intervals were checked for missing data, which may corrupt the time series. Sporadic and random missing data do not appear in groups and therefore do not corrupt the time series after interpolation. If any of those intervals had more than 50% of data missing, the interval was excluded for all measurement levels. The remaining intervals were tested for Taylor’s hypothesis (TH), which allows us to transform from the space domain (wavenumber, k) to the time domain (frequency, f ). The criterion for the validity of TH is σ < 0.5U, where U is the mean horizontal wind speed and σ is the standard deviation (e.g., [44]). If the ratio of the standard deviation to mean horizontal wind speed was greater than or equal to 0.5, the corresponding intervals were omitted.


7 of 25

2.6. Stability Parameter and Friction Velocity In order to analyze thermal stratification of each bora episode, dimensionless height ζ, or stability parameter (e.g., [44]) was calculated for each block interval using the following equation: ζ=

z L

(1)

Here, z represents the height of the observation level and L is the Obukhov length defined as the height where dynamical-mechanical turbulence generation is approximately equal to the thermal-buoyancy contribution to the turbulence generation or destruction: L=

−θu∗ 3 gk(w0 θ 0 )

(2)

where k = 0.4 is the von Karman constant (e.g., [47]), g = 9.81 m·s−2 is acceleration due to gravity, θ is the mean virtual potential temperature, w0 θ 0 is local (at each measurement height) vertical kinematic heat flux, and u∗ is the local friction velocity calculated as: q u∗ =

4

2

(u0 w0 ) + (v0 w0 )

2

(3)

Since the relative humidity measurements were not available, we used ultrasonic temperature values instead of potential virtual temperature for determining turbulence heat flux. According to some authors [30,48,49], the ultrasonic temperature is a good approximation of the potential temperature. The negative stability parameter implies statically unstable stratification; a positive stability parameter implies stable stratification, while the stratification is neutral when the stability parameter is equal to zero. Here, we took block intervals with absolute values of the stability parameter less than 0.02 as near-neutral, which is the same criterion that was used in [30]. 2.7. Monin–Obukhov Similarity Functions Due to the intensive mechanical mixing during bora episodes, the thermal stratification of the atmosphere is generally very close to neutral (ζ ≈ 0) [35]. However, there are periods when the value of the stability parameter deviates from zero, which means that in those periods the atmosphere is not neutrally stratified. According to the Monin–Obukhov similarity theory, during those periods, the wind shear should satisfy the following relationship ∂U u∗ = Φm (ζ ) ∂z kz

(4)

where U is the mean value of the streamwise wind speed (u = U + u0) and Φm (ζ ) represents the following similarity function: Φm (ζ ) = 1 + 4.7ζ for ζ > 0 (5) Φm (ζ ) = (1 − 15ζ )−1/4

for ζ < 0

(6)

In order to check the applicability of the similarity theory in bora wind cases, the block intervals were divided according to the thermal stratification of the atmosphere into stable and unstable intervals. The experimental value of the similarity function appropriate for each block interval was calculated using the measurement data as follows: Φm, exp (ζ 2 ) =

∂U kz ∂z u∗2

(7)

where Φm, exp is the experimental value of the similarity function, ζ 2 is the stability parameter at the 2 m height, and u∗2 is the friction velocity, also at the 2 m height. Since the measurement data were


8 of 25

available at measurement heights of 2, 5, and 10 m, two layers were analyzed: 2–5 m and 5–10 m. The height (z) was calculated as the arithmetic mean for each layer and the wind shear (∂U/∂z) was calculated with finite differences. 2.8. Turbulence Kinetic Energy Budget TKE is a measure of turbulence intensity. Therefore, it is one of the most important variables in micrometeorology. The equations needed to examine the TKE budget and individual terms are implemented following Equations (1) and (2) in [49]. The individual terms in the TKE budget equation (e.g., [44]) describe the physical processes that generate, transport, or suppress turbulence. After rotating the coordinate system in the mean wind direction, V is globally zero for each episode. The standard assumption is that there is no subsidence (W = 0). We checked and confirmed that W was small compared to U and could be neglected. The one-dimensional, nearly horizontally homogeneous TKE budget equation can be written as (hereafter in the text e¯ will be referred to as TKE for simplicity): ∂ w0 e g ∂e ∂U = + w0 θ 0 − u0 w0 − − ε + Rs (8) ∂t ∂z ∂z θ I

II

III

IV

where U is the mean value of the streamwise wind speed, u’ and w’ are the corresponding turbulent values. Variable ε represents the dissipation of TKE into heat by molecular viscosity, and Rs is the residual. In theory, the overbars represent suitable spatial (instead of ensemble) averaging. Since we require the validity of TH, in practice, the overbars are considered as time averaging on block intervals. Term I represents a local storage of TKE, which is equal to the increase or decrease of TKE in time at a given location due to all of the TKE production, transport, or redistribution and destruction terms. The production and destruction terms include the buoyant production/consumption (term II), which depends on the sign of the heat flux; the mechanical (shear) production (term III), which is typically positive in the surface layer because of opposite signs of the horizontal momentum fluxes and vertical wind shear; the vertical turbulent transport or redistribution of TKE by turbulent eddies (term IV); the dissipation of TKE into heat by molecular viscosity (ε); and the residual Rs containing term that describes the pressure transport of TKE. All the TKE terms, besides the dissipation rate (ε), can easily be calculated directly from the measured u, v, and w components. The local change of TKE is calculated using the central finite difference scheme. Term II is calculated for 2 m, 5 m, and 10 m separately, and the middle level (3.5 m and 7.5 m) values are obtained by averaging between the two heights. Similarly, all terms that require values from upper and lower levels are calculated for the corresponding middle level using the central finite difference scheme. In order to evaluate ε, the inertial dissipation method (IDM)—provided by Kolmogorov’s 1941 hypothesis—can be employed following Equations (7)–(9) in [50]. The same method was used by Veˇcenaj in [33,34]. Finally, the residual term is calculated using all known terms by assuming the balance between the left- and right-hand sides of Equation (8). 2.9. Wind Profiles In order to investigate the agreement of the experimental data with the logarithmic-law approximation for neutral wind speed vertical profiles in the surface layer [44], statically near-neutral block intervals during bora episodes were studied. These intervals are defined as the intervals during which |ζ | < 0.02 is valid for all three levels [30]. For every interval the profile friction velocity u∗p and roughness length z0 are calculated as follows: u ∗p =

k (U10 − U2 ) ln zz102

(9)


9 of 25

z0 = z2 exp −

U2 z ln 10 U10 − U2 z2

(10)

where z10 and z2 are the heights, while U10 and U2 are time-averaged mean wind speeds in the x-direction at the highest (10 m) and the lowest (2 m) levels. The vertical wind profile is reconstructed with mean and median values of these parameters as: u(z) =

u ∗p z ln k z0

(11)

Equations (9)–(11) are derived from the Equation (3), using certain approximations and parameterizations [51]. Thus, they are related to the mean, low-frequency measurements, in contrast to local friction velocity, which is related to high-frequency measurements. 3. Results and Discussion 3.1. Selected Bora Episodes Table 2 shows all the detected bora episodes, classified by their types. Synoptic maps and soundings analyses indicate that four episodes are transitional with respect to flow depth (B02, B07, B08, and B11 in Table 2), while in two episodes (B03 and B12) the synoptic flow changes from cyclonic to anticyclonic without significant change of the flow depth. Three episodes (B05, B10, and B14) also include passages of occluded or warm fronts, but none of the episodes exhibit typical characteristics of frontal bora type. Episode B05 is very short and partially caused by cold air advection as the surface cyclone progressed southeast, but the surface analysis did not show the corresponding surface cold front. The shortest and also the weakest episode (B10) has some characteristics of frontal bora, although surface analysis did not show a typical cold front passage, but merely a slow transition of a quasi-stationary front, as the cold air mass slowly advected from the northwest. Due to the low wind speeds and untypical character of these episodes, they were not included in our micro-scale turbulence analysis. Table 2. All detected bora episodes. Times are in UTC. “Avg U” is the average 10-min wind speed at 10 m, “max U” is the maximum 10-min wind speed value, and “max G” is the 1-s maximum gust. All wind speeds are in m·s−1 . The chosen episodes and dominant bora type are marked with bold font. Episode

Start

End

Duration

Avg U

Max U

Max G

Type

B01 B02 B03 B04 B05 B06 B07 B08 B09 B10 B11 B12 B13 B14

10 October 2015. 09:20 21 October 2015. 16:00 23 October 2015. 11:30 30 October 2015. 03:30 21 November 2015. 15:50 22 November 2015. 08:20 26 November 2015. 01:30 10 Decemebr 2015. 01:40 30 Decemebr 2015. 10:10 3 January 2016. 12:40 15 January 2016. 23:30 17 January 2016. 14:50 21 January 2016. 13:00 3 February 2016. 22:40

11 October 2015. 07:10 23 October 2015. 01:40 24 October 2015. 00:30 1 November 2015. 12:30 22 November 2015. 00:40 22 November 2015. 20:10 28 November 2015. 04:00 10 Decemebr 2015. 12:00 30 Decemebr 2015. 21:40 3 January 2016. 19:00 17 January 2016. 04:10 18 January 2016. 07:30 22 January 2016. 06:40 4 February 2016. 16:50

0 days 21:50 1 day 09:40 0 days 13:00 2 days 09:00 0 days 08:50 0 days 11:50 2 days 02:30 0 days 10:20 0 days 11:30 0 days 06:20 1 day 04:40 0 days 16:40 0 days 17:40 0 days 18:10

12.16 13.35 8.08 7.16 10.31 7.44 17.69 16.57 9.58 5.92 18.02 20.82 7.25 19.17

18.21 20.12 13.53 12.88 16.77 12.26 31.05 25.66 16.46 9.62 28.40 28.78 10.91 25.83

27.94 35.68 24.18 21.02 26.34 18.60 45.47 36.10 26.90 16.35 40.18 39.68 17.42 37.27

SC SC/DC DC/DA SA F/SC SC SC/DC SA/DA DA SC/F SC/DC SC/SA SA F/SC

Bora B14 was predominantly of a shallow cyclonic (SC) type with a warm front passage and cold air advection following the cyclone. Non-transitional deep cyclonic (DC) bora was not detected in the recorded data. Thus, for further micro-scale analysis, we chose the following episodes: B01 (SC), B09 (DA), and B13 (SA), because they represent typical cases throughout the episode. Although B04 (SA) was longer and had a higher maximum wind speed than B13 (SA), the sounding data were missing for part of the episode, so B13 (SA) was chosen instead. Episode B13 actually had a higher


10 of 25

average wind speed than B04. Additionally, note from Table 2 that 1-s gust maxima are usually two to three times larger than the related 10-min speed averages at 10 m. Figure 3 shows the surface analysis (NCEP GDAS/FNL) at 12 UTC on 10 October 2015, 2018, 9, PEER 10 Atmosphere 2018, the 9, xx FOR FOR PEER REVIEW REVIEW 10 of of 25 25 whichAtmosphere represents synoptic situation after the beginning of the B01 (SC) episode.

Figure 3. The surface analysis (NCEP(NCEP GDAS/FNL) mean sea level pressure (solid line) andline) geopotential Figure Figure 3. 3. The The surface surface analysis analysis (NCEP GDAS/FNL) GDAS/FNL) mean mean sea sea level level pressure pressure (solid (solid line) and and height at 500 hPa (dashed line) for 12 UTC on 10 October 2015. The bora episode B01 (SC). geopotential geopotential height height at at 500 500 hPa hPa (dashed (dashed line) line) for for 12 12 UTC UTC on on 10 10 October October 2015. 2015. The The bora bora episode episode B01 B01 The measurement site is marked with the orange star. star. (SC). site with the (SC). The The measurement measurement site is is marked marked with the orange orange star.

The surface cyclone, with itsits theTyrrhenian TyrrhenianSea, Sea, influences eastern The cyclone, with center situated over influences the eastern The surface surface cyclone, with itscenter centersituated situated over over the the Tyrrhenian Sea, influences the the eastern Adriatic coast coast with its northeastern quadrant. A more analysis of the surface and upper-level Adriatic with its quadrant. A more analysis of and Adriatic coast with its northeastern northeastern quadrant. A detailed more detailed detailed analysis of the the surface surface and upper-level features reveal strong surface gradients Alps an features reveal strong surface over the Dinaric over Alps the andDinaric an upper-level trough upper-level features revealpressure strong gradients surface pressure pressure gradients over the Dinaric Alps and and an at upper-level trough at with winds the This is typical upper-level trough at 500 500 hPa withThis winds from the south. south. This for is aa SC typical situation for SC SC bora bora type. type. 500 hPa with winds from thehPa south. is from a typical situation borasituation type. for The 12 sounding for the same day ideally upstream from the The Zagreb Zagreb 12 UTC UTC soundingfor forthe thesame same day day (Figure (Figure 4a) isisnot not ideally upstream fromfrom the the The Zagreb 12 UTC sounding (Figure4a) 4a)is not ideally upstream measuring tower (it is more to the north compared to low-level wind, which is more from the east), measuring tower (it is more to the north compared to low-level wind, which is more from the east), measuring tower (it is more to the north compared to low-level wind, which is more from the east), but it still the structure of bora conditions in lower levels. At the same still reflects reflects the vertical vertical structure of upstream upstream bora conditions lowerlevels. levels.At Atthe thesame sametime, but itbut stillitreflects the vertical structure of upstream bora conditions inin lower time, the Zadar sounding (Figure 4b) shows that the bora layer ends above 850 hPa, which time, the Zadar sounding (Figure 4b) shows that the bora layer ends above 850 hPa, which the Zadar soundingan(Figure 4b) showsatthat the bora layer ends above 850 hPa, which corresponds to corresponds corresponds to to an inversion inversion visible visible at around around 800 800 hPa. hPa. an inversion visible at around 800 hPa.

(a) (a)

(b) (b)

Figure 4. skew-T log-P graph of sounding the Zagreb (a) (b) at Figure 4. The The skew-T log-Pgraph graphof of sounding sounding data data from thethe Zagreb (a) and and Zadar (b) stations stations at 12 12 at Figure 4. The skew-T log-P datafrom from Zagreb (a) Zadar and Zadar (b) stations UTC on 10 October 2015. Bora episode B01 (SC). The vertical axis (pressure) is in hPa and horizontal UTC on 10 October 2015. Bora episode B01 (SC). The vertical axis (pressure) is in hPa and horizontal 12 UTC on 10 October 2015. Bora episode B01 (SC). The vertical axis (pressure) is in hPa and horizontal axis is in °C. axis (temperature) (temperature) axis (temperature) is inis◦in C.°C.

Atmosphere Atmosphere 2018, 2018, 9, 9, 116 x FOR PEER REVIEW Atmosphere 2018, 9, x FOR PEER REVIEW

11 of 25 11 of 25

Figure 5a shows the time series of 10 Hz streamwise wind speed (u) at the 10 m level for the Figure 5a 5a shows shows the the time time series series of of 10 10 Hz Hz streamwise streamwise wind wind speed speed (u) (u) at at the the 10 10 m m level level for for the the Figure bora B01. bora B01. bora B01.

Figure 5. (a) Bora episode (B01) streamwise wind speed (u) downsampled to 10 Hz. The white line isis Figure5. 5.(a) (a)Bora Boraepisode episode(B01) (B01)streamwise streamwisewind windspeed speed(u) (u)downsampled downsampledto to10 10Hz. Hz.The Thewhite whiteline lineis Figure the 10-min moving average. (b) The zoomed 60 min of the same data with 1-min moving average is the10-min 10-minmoving movingaverage. average.(b) (b)The Thezoomed zoomed6060 min same data with 1-min moving average the min ofof thethe same data with 1-min moving average is inis in color. in red red color. red color.

The zoomed wind speed data with aa 1-min moving average (Figure 5b) the Thezoomed zoomedwind windspeed speed data with 1-min moving average (Figure 5b) shows shows the pulsations pulsations The data with a 1-min moving average (Figure 5b) shows the pulsations with −1 and varying periods on a scale of around 5 min. with an amplitude of more than 10 m· s − 1 −1 and varying periods on a scale of around 5 min. with an amplitude of more than 10 m· s an amplitude of more than 10 m·s and varying periods on a scale of around 5 min. For episode B09 (DA) surface analysis at 12 on 30 2015 6) Forepisode episodeB09 B09(DA) (DA) surface analysis at UTC 12 UTC UTC onDecember 30 December December 2015 (Figure (Figure 6) features features For surface analysis at 12 on 30 2015 (Figure 6) features a strongaa strong anticyclone (1045 hPa at the center) influencing most of Europe (a larger-area surface analysis strong anticyclone (1045 hPacenter) at the influencing center) influencing of Europe (a larger-area anticyclone (1045 hPa at the most ofmost Europe (a larger-area surface surface analysisanalysis can be can be found at [40]). can be found at [40]). found at [40]).

30 December 2015. Figure Figure6.6. As As with withFigure Figure3,3,but butfor forB09 B09(DA); (DA);12 12UTC UTCon on30 30December December2015. 2015.

At the same time, which is approximately 2 h after the beginning of the B09 episode, 500 hPa At the same same time, time, which is approximately 2 h after the beginning beginning of the B09 episode, episode, 500 hPa geopotential height (Figure 6) shows the western flank of an upper-level trough, passing over the upper-level trough, passing geopotential height (Figure 6) shows the western flank of an upper-level trough, passing over the Black Black Sea Sea region. region. This This combination combination of of the the surface surface and and upper upper level level features features provides provides N/NE N/NE winds winds throughout the troposphere. throughout the troposphere.


12 of 25

Black Sea2018, region. This combination Atmosphere 9, x FOR PEER REVIEW Atmosphere 2018, 9, x FOR PEER REVIEW

of the surface and upper level features provides N/NE winds 12 of 25 12 of 25 throughout the troposphere. Deep flow can can also also be seen in the Zagreb sounding Deep N/NE N/NE flow sounding (Figure 7a), 7a), together together with with strong strong Deep N/NE flow can also be seen in the Zagreb sounding (Figure 7a), together with strong subsidence 800800 hPa—a typical signature of such stronga anticyclone. The ZadarThe sounding subsidenceinversion inversionat at hPa—a typical signature ofa such strong anticyclone. Zadar subsidence inversion at 800 hPa—a typical signature of such a strong anticyclone. The Zadar (Figure 7b) features a similar wind profile, with somewhat weaker subsidence inversion. sounding (Figure 7b) features a similar wind profile, with somewhat weaker subsidence inversion. sounding (Figure 7b) features a similar wind profile, with somewhat weaker subsidence inversion.

(a) (a)

(b) (b)

Figure 7. As with Figure 4, but for the bora episode B09 (DA); 12 UTC on 30 December 2015. Figure As with Figure 4,for butthe forbora the bora episode B09 (DA); 12 UTC 30 December Figure 7. As7.with Figure 4, but episode B09 (DA); 12 UTC on 30onDecember 2015.2015.

Figure 88 shows shows the the time time series of of streamwisewind wind speed(u) (u) atthe the 10m m levelfor for boraB09. B09. Figure Figure 8 shows the time series series of streamwise streamwise wind speed speed (u) at at the 10 10 m level level for bora bora B09.

Figure 8. As with Figure 5, but for the bora episode B09 (DA). Figure Figure 8. 8. As As with with Figure Figure 5, 5, but but for for the the bora bora episode episode B09 B09 (DA). (DA).

The bora episode B09 has somewhat lower wind speeds in the middle and higher wind speeds The The bora bora episode episode B09 B09 has has somewhat somewhat lower lower wind wind speeds speeds in in the the middle middle and and higher higher wind wind speeds speeds after the beginning and before the weaker end of the episode. The zoomed data (Figure 8b) with after beginning and the weaker weaker end end of of the the episode. episode. The after the the beginning and before before the The zoomed zoomed data data (Figure (Figure 8b) 8b) with with 1-min moving average shows the superposition of pulsations with varying periods (around 3–5 1-min moving average shows the superposition of pulsations with varying periods (around 3–5 1-min moving average shows the superposition of pulsations with varying periods (around 3–5 min). min). Some of the pulsations have estimated amplitudes over 10 m·s−1−1. − 1 min). the pulsations have estimated amplitudes over Some Some of the of pulsations have estimated amplitudes over 10 m·s10 m· . s . The surface analysis at 00 UTC on 22 January 2016 (Figure 9), 11 h after the beginning of the B13 The surface analysis at 00 UTC on 22 January 2016 (Figure 9), The surface analysis at 00 UTC on 22 January 2016 (Figure 9), 11 11 h h after after the the beginning beginning of of the the B13 B13 episode (SA), shows a weak anticyclone with the center over Central Europe. episode (SA), shows a weak anticyclone with the center over Central Europe. episode (SA), shows a weak anticyclone with the center over Central Europe.

Atmosphere 2018, 9, 116 Atmosphere 2018, 9, x FOR PEER REVIEW Atmosphere 2018, 9, x FOR PEER REVIEW

13 of 25 13 of 25 13 of 25

Figure 9. 9. As with Figure 3, but but for for B13 B13 (SA); (SA); 00 00 UTC on on 22 January January 2016. Figure As with with Figure Figure 3, 3, Figure 9. As but for B13 (SA); 00 UTC UTC on 22 22 January 2016. 2016.

moredetailed detailedanalysis analysis shows that the strongest surface pressure gradients are over the A more shows thatthat the the strongest surface pressure gradients are over theover Dinaric A more detailed analysis shows strongest surface pressure gradients are the Dinaric Alps, while the 500 hPa NW (northwesterly) winds at the western flank of the upper-level Alps, while the 500 hPa NW (northwesterly) winds at the western flank of the upper-level trough Dinaric Alps, while the 500 hPa NW (northwesterly) winds at the western flank of the upper-level trough dominate dominate the area. area. dominate the area.the trough January 2016 The Zagreb at 00 00 UTC The Zagreb sounding sounding at UTC on on 22 22 January January 2016 (Figure (Figure 10a) 10a) displays displays aa shallow shallow NE NE wind wind with NW winds throughout the rest of the troposphere. layer, capped by an inversion at 850 hPa, throughout the rest of the troposphere. layer, capped by an inversion at 850 hPa, with NW winds throughout the rest of the troposphere. Zadar sounding sounding (Figure (Figure 10b) 10b) also also shows shows aa shallow shallow NE NE wind wind layer, layer, capped capped by by aa very very stable stable layer layer and and Zadar with NW NW winds winds above. above. with with NW winds above.

(a) (a)

(b) (b)

Figure 10. As with Figure 4, but for B13 (SA); 00 UTC on 22 January 2016. Figure 10. As with Figure 4, but for B13 (SA); 00 UTC on 22 January 2016.

The time time series series of streamwise wind speed (u) atatthe the 1010m mmlevel level forfor this episode can be be seen in series of of streamwise streamwisewind windspeed speed(u) (u)at the10 levelfor this episode can seen The this episode can be seen in Figure 11. in Figure Figure 11.11.


14 of 25

Atmosphere 2018, 9, 9, x FOR PEER REVIEW Atmosphere 2018, x FOR PEER REVIEW

1414 ofof 2525

Figure 11. AsAs with Figure 5,5, but for the bora episode B13 (SA). Figure 11. As with Figure 5, but for the bora episode B13 (SA). Figure 11. with Figure but for the bora episode B13 (SA).

As Asinin inthe theprevious previoustwo twoepisodes, episodes,a aa1-min 1-minmoving movingaverage averageofof ofthe thezoomed zoomeddata data(Figure (Figure11b) 11b)also also As the previous two episodes, 1-min moving average the zoomed data (Figure 11b) also 1 and shows s·−1 periods around 5 5min. showspulsating pulsatingbehavior behaviorwith withestimated estimatedamplitudes amplitudesaround around555m· m· periods around shows pulsating behavior with estimated amplitudes around m ss−−1and and periods around 5 min. min. 3.2. Bora Spectra 3.2. 3.2.Bora BoraSpectra Spectra Frequency weighted power spectral densities were calculated and smoothed (Section 2.4) for Frequency Frequencyweighted weightedpower powerspectral spectraldensities densitieswere werecalculated calculatedand andsmoothed smoothed(Section (Section2.4) 2.4)for for horizontal wind components at levels 2, 5, and 10 m for all episodes. In Figure 12, the spectra are horizontal horizontalwind windcomponents componentsat levels 2, 5,5,and and1010mmfor forallallepisodes. episodes.In Figure 12, the thespectra spectraare are shown for three selected bora episodes (B01, B09, and B13). shown shownfor forthree threeselected selectedbora boraepisodes episodes(B01, (B01,B09, B09,and andB13). B13).

Figure 12. log-linearrepresentation representationofofthe the frequency-weighted power spectral densities the Figure Figure12. 12.AAlog-linear thefrequency-weighted frequency-weightedpower powerspectral spectraldensities densitiesofof ofthe the longitudinal (u) component (a,c,e) and lateral (v) component (b,d,f) of the wind speed at levels 2, longitudinal longitudinal(u)(u)component component(a,c,e) and lateral (v) component (b,d,f) of the wind speed at levels 2, 5,5, and m, for the three selected bora episodes; (a,b) B01 (SC); (c,d) B09 (DA); (e,f) B13 (SA). Thick and and1010 10m, m,for forthe thethree threeselected selectedbora boraepisodes; episodes;(a,b) (a,b)B01 B01(SC); (SC);(c,d) (c,d)B09 B09(DA); (DA);(e,f) (e,f)B13 B13(SA). (SA).Thick Thick dashed black vertical lines indicate the 30-min and 5-min periods. dashed black vertical lines indicate the 30-min and 5-min periods. dashed black vertical lines indicate the 30-min and 5-min periods.

clearly visible uucomponent contains more energy than thethe vthe component for all three ItIt that the contains more energy than v vcomponent for allall Itisisisclearly clearlyvisible visiblethat thatthe the ucomponent component contains more energy than component for types of bora (note the y-axis ranges), which was expected due to the initial rotation of the coordinate three threetypes typesofofbora bora(note (notethe they-axis y-axisranges), ranges),which whichwas wasexpected expecteddue duetotothe theinitial initialrotation rotationofofthe the coordinate coordinatesystem. system.Additionally, Additionally,asasexpected, expected,the thedata dataatat1010mmcontain containa ahigher higheramount amountofofenergy energy than thanthe thedata dataatat5 5and and2 2mmatatthe thelower lowerfrequency frequencyband. band.However, However,the thedata dataatat2 2mmcontains containsa ahigher higher


15 of 25

Atmosphere 2018, 9, x FOR PEER REVIEW

15 of 25

system. Additionally, as expected, the data at 10 m contain a higher amount of energy than the data at 5amount and 2 mofatenergy the lower band. However, the the datadata at 2 at m contains higher amountin of spectra energy at afrequency higher frequency band than 10 and 5am. Differences at a higher band than the at 10 m. Differences in spectra between u and between thefrequency u and v components aredata larger forand the5 bora episodes with stronger winds the (B01 and vB09), components are larger for the bora episodes with stronger winds (B01 and B09), depending on a depending on a frequency band in which the energy is observed. frequency band in which the energy is observed. The appropriate turbulence averaging time scale had to be estimated in order to proceed with The appropriate turbulence averaging time chosen scale had to be estimated in order proceed with the analysis (Section 2.4). In this study, we have to do it in a subjective way,to looking only at the analysis (Section 2.4). In this study, we have chosen to do it in a subjective way, looking only at the power spectral densities and bearing in mind some most recent bora turbulence research [1,30]. the power spectral densities and“gap”) bearing mind some most recent turbulence [1,30]. A minimum of energy (energy in in frequency-weighted powerbora spectra for boraresearch episodes B01 A minimum of energy (energy “gap”) in frequency-weighted power spectra for bora episodes B01 (SC) and B09 (DA) is located near the frequency corresponding to the 30-min period. Spectra for (SC) B09 (DA) is located near the afrequency corresponding to did the 30-min for bora and episode B13 (SA) do not show clear minimum like they for B01period. and B09Spectra or, more bora episode (SA)from do not show a clear minimum likeside theyofdid B01 and is B09 or, more precisely, precisely, theB13 energy large-scale motions (the left thefor spectrum) missing. Regardless the energy from large-scale motions (the left side of the spectrum) is missing. Regardless of the of the missing energy, there is a significant difference in the energy between periods of 30missing and 5 energy, there is a significant difference in the energy periods 30 and 5 min,B13. and we can say min, and we can say that there is also an energy gapbetween in the spectra forofbora episode that there alsoenergy an energy in the spectra for15 bora B13. Sinceisthe gap gap is located between andepisode 40 min, we may use a 30-min time period Since the energy gap is located between 15 and 40 min, we may with use a[52] 30-min time period following the instructions provided in Section 2.5. This is in accordance for the nocturnal following the instructions provided in Section Thisofis Kutina), in accordance for [30] the nocturnal stable boundary layer in the Croatian lowland2.5. (town while with Babić[52] et al. found an stable boundary layer in the Croatian lowland (town of Kutina), while Babi´ c et al. [30] found an energy gap at the 15-min period for bora episodes at Pometeno Brdo, northeast from the energy city of gap the 15-min periodanalyses for borawere episodes at Pometeno Brdo, northeast from city of Split. Hence, Split.atHence, all further performed on block intervals of 30 minthe length. all further werepeak performed block intervals of 30 min Thereanalyses is an energy in bothonthe u and v component of length. the given power spectral densities, There is an energy peak in both the u and v component of the given powerrelated spectraltodensities, between between the periods of ≈2 and 8 min. The energy peaks are most likely bora pulsations the periodswhich of ≈2 can and also 8 min. most(Figures likely related to bora [11,28,29], be The seenenergy in the peaks u timeare series 5, 8, and 11). pulsations The power[11,28,29], spectral which seen inB09 the u time series (Figures 5, 8 and periods 11). The 104 power spectral densitycan of also borabeepisode shows a large peak between s and 103 s.density This isofatbora the episode B09 shows a large peak between periods 104 s and 103 s. This is at the lowest end of the lowest end of the pulsations similar to that in [53]. Spectral analysis did not show any significant pulsations similar that in [53]. Spectral analysis not show significant difference power difference in powertospectral densities between boradid episodes B01any (SC) and B13 (SA). Otherin than the spectral densities between bora episodes B01 (SC) and B13 (SA). Other than the peak in the B09 (DA) peak in the B09 (DA) spectra, which could be related to different flow dynamics (deep bora), there spectra, which could be related the to different flow dynamics bora), there is no clearand connection is no clear connection between power spectral densities(deep of different bora episodes the type between of bora. the power spectral densities of different bora episodes and the type of bora. 3.3. Friction Velocity Velocity and 3.3. Friction and Stability Stability Parameter Parameter Time series of friction velocity Time series of friction velocity calculated calculated from from turbulent turbulent fluxes fluxes (Section (Section 2.6, 2.6, Equation Equation (3)), (3)), for for all all three The friction friction velocity velocity time time series series is is closely closely related related three levels levels of of bora bora B01 B01 (SC) (SC) are are shown shown in in Figure Figure 13. 13. The to to the the TKE TKE time time series series (not (not shown). shown).

Figure 13. The time series of friction velocity (𝑢∗ ) for bora B01 (SC). The blue line is 𝑢∗ at 2 m, the Figure 13. The time series of friction velocity (u∗ ) for bora B01 (SC). The blue line is u∗ at 2 m, the green green line is at 5 m, and the magenta line is at 10 m. line is at 5 m, and the magenta line is at 10 m.

Figure 13 shows that lower values of friction velocity appear at the beginning and at the end of the bora episode. Comparing this figure to Figure 5, it can be seen that the friction velocity is closely related to the streamwise wind speed (u). It should also be noted that the friction velocity values are


16 of 25

Figure 13 shows that lower values of friction velocity appear at the beginning and at the end of the bora episode. Comparing this figure to Figure 5, it can be seen that the friction velocity is closely Atmosphere 2018, 9, x FOR PEER REVIEW 16 of 25 Atmosphere x FOR PEER wind REVIEW 16 ofare 25 related to2018, the 9, streamwise speed (u). It should also be noted that the friction velocity values generally highest at the 5 m level. In total, the bora B01 (SC) has 43 30-min blocks, with all blocks generally highest at the 5 m level. In total, the bora B01 (SC) has 43 30-min blocks, with all blocks generally highest at the 5 m level. In total,2.5). the bora B01 (SC) has 43 30-min blocks, with all blocks satisfying satisfying Taylor’s Taylor’s hypothesis hypothesis (TH, (TH, Section Section 2.5). satisfying Taylor’s hypothesis (TH, Section 2.5). The B09 The B09 bora bora (DA) (DA) (Figure (Figure 14) 14) is is relatively relatively short, short, and and thus thus has has only only 23 23 30-min 30-min blocks, blocks, of of which which The B09atbora (DA)of(Figure 14) isdid relatively short, and thus has only 23invalid 30-minblocks blocks, ofusually which three blocks the end the episode not satisfy TH (Section 2.5). Such are three blocks at the end of the episode did not satisfy TH (Section 2.5). Such invalid blocks are usually three blocks at the end of the episode did not satisfy TH (Section 2.5). Such invalid blocks are usually located the beginning beginning or the end end of of the the episode, episode, where where turbulence turbulence intensity intensity seems be the the highest highest located at at the or the seems to to be located at standard the beginning or the end of the low episode, where turbulence intensity seems to bethat the highest (i.e., large deviation of relatively wind speed). Nevertheless, Figure 14 shows friction (i.e., large standard deviation of relatively low wind speed). Nevertheless, Figure 14 shows that (i.e., large standard deviation of relatively lowthis wind speed). Nevertheless, Figure 14 explains shows that velocity is closely related to the wind speed for episode as well (Figure 8). This also the friction velocity is closely related to the wind speed for this episode as well (Figure 8). This also friction velocity is closely related to the wind speed for this episode as well (Figure 8). This also lower friction velocity values in the values middleinofthe thismiddle episode, as the wind speed also lower in also that explains the lower friction velocity of this episode, as thewas wind speed was explains the lower friction velocity values in the middle of this episode, as the wind speed was also part of the episode. lower in that part of the episode. lower in that part of the episode.

Figure 14. As with Figure 13, but for the bora episode B09 (DA). Figure 14. As with Figure 13, but for the bora episode B09 (DA).

Interestingly, forthis this episode,the the friction velocity at the m level is also higher than the Interestingly, for friction velocity at the 5 m5 is also higher than at theat other Interestingly, for thisepisode, episode, the friction velocity at the 5level m level is also higher than at the other two levels. two otherlevels. two levels. Bora B13 (SA), 15, has has 35 of which one block the end Bora B13 (SA), in in Figure Figure 15, 35 30-min 30-min blocks, blocks, of which one block near near the end of of the the episode episode Bora B13 (SA), in Figure 15, has 35 30-min blocks, of which one block near the end of the episode was discarded. was discarded. was discarded.

Figure 15. As with Figure 13, but for the bora episode B13 (SA). Figure 15. As with Figure 13, but for the bora episode B13 (SA). Figure 15. As with Figure 13, but for the bora episode B13 (SA).

As with the other two episodes (B01 and B09), a close relation between wind speed (Figure 11) As with the other two episodes (B01 and B09), a close relation between wind speed (Figure 11) and friction velocity can easily be seen. For this episode, the friction velocity at 5 m is also higher and friction velocity can easily be seen. For this episode, the friction velocity at 5 m is also higher than at the other two levels. than at the other two levels. The distribution of friction velocity for each bora episode at different vertical levels is shown in The distribution of friction velocity for each bora episode at different vertical levels is shown in Figure 16. For the bora episode B01 (SC), two maxima occur at the lower levels. The first one (located Figure 16. For the bora episode B01 (SC), two maxima occur at the lower levels. The first one (located around 0.7–0.8 m·s−1 ) probably corresponds to the moderate wind speeds at the beginning and the −1 around 0.7–0.8 m·s ) probably corresponds to the moderate wind speeds at the beginning and the


17 of 25

As with the other two episodes (B01 and B09), a close relation between wind speed (Figure 11) and friction velocity can easily be seen. For this episode, the friction velocity at 5 m is also higher than at the other two levels. The distribution of friction velocity for each bora episode at different vertical levels is shown in Atmosphere x FOR PEER REVIEW 17 of 25 Figure 16.2018, For9,the bora episode B01 (SC), two maxima occur at the lower levels. The first one (located − 1 around 0.7–0.8 m·s ) probably corresponds to the moderate wind speeds at the beginning and the end of the episode, while the other one, which is larger (located around 0.9–1 m·s−−11), could be a result end of the episode, while the other one, which is larger (located around 0.9–1 m·s ), could be a result of the wind speed increase during the more developed part of the episode. of the wind speed increase during the more developed part of the episode.

Figure 16. Frequency distribution of friction velocity for: the B01 (a–c), B09 (d–f), and B13 (g–i) bora Figure 16. Frequency distribution of friction velocity for: the B01 (a–c), B09 (d–f), and B13 (g–i) bora episodes; at three different vertical levels: 2 m (a,d,g), 5 m (b,e,h), and 10 m (c,f,i). episodes; at three different vertical levels: 2 m (a,d,g), 5 m (b,e,h), and 10 m (c,f,i).

This bimodality is less pronounced at the highest level, where the friction velocity distribution This bimodality is less pronounced at the highest level, where the friction velocity distribution tends to look more near-normal with similar mean and median values (Table 3). This shift towards tends to look more near-normal with similar mean and median values (Table 3). This shift towards near-normal distribution is in better agreement with the results found in [35]. near-normal distribution is in better agreement with the results found in [35]. Other3.two bora of episodes seem to have such shifted towards The anticyclonic summary statistics the friction velocity u∗ (m ·s–1 ) distribution for each bora types episode at selected Table lowervertical valueslevels. of friction velocity (maximum located around 0.5 m·s−1) compared to bora B01 (SC). A slight negative skewness can be noticed for those bora episodes, and is also evident as somewhat larger medians compared to the mean velocity. Min Episode Height (m)corresponding Mean MedfrictionMax SD 2 0.782 0.830 1.051 0.384 0.177 Table 3.B01 The summary statistics of the friction velocity s–1) for each bora episode at selected 5 0.838 0.832 𝑢∗ (m· 1.118 0.465 0.171 (SC) vertical levels. 10 0.739 0.737 1.051 0.299 0.186 0.592Med 0.892 0.247 SD 0.177 Height 0.580 (m) Mean Max Min 0.656 0.681 0.999 0.246 0.214 2 0.435 0.7820.4210.830 0.679 1.051 0.384 0.1710.177 0.154 B01 (SC) 5 0.838 0.832 1.118 0.465 0.171 2 0.454 0.477 0.673 0.106 0.138 10 0.497 0.7390.5300.737 0.701 1.051 0.299 5 0.1530.186 0.135 B13 (SA) 2 0.419 0.5800.4510.592 0.622 0.892 0.247 10 0.1710.177 0.120 B09 (DA) 5 0.656 0.681 0.999 0.246 0.214 10 0.435 0.421 0.679 0.171types 0.154 Other two anticyclonic bora episodes seem to have such distribution shifted towards lower 2 0.454 0.477 0.673 0.106 0.138 − 1 values of friction velocity (maximum located around 0.5 m·s ) compared to bora B01 (SC). A slight 5 those bora 0.497episodes, 0.530 and 0.701is also 0.153 0.135as somewhat larger negative skewness B13 can (SA) be noticed for evident 10 0.419 0.451 0.622 0.171 0.120 medians compared to the corresponding mean friction velocity. B09 (DA)

Episode2

5 10

The stability parameter distribution is shown in Figure 17. In order to compare the stability The stability parameter distribution is shown in Figure 17. In order to compare the stability parameters between different bora episodes, the width of the bin was kept the same (0.02). parameters between different bora episodes, the width of the bin was kept the same (0.02).

Atmosphere Atmosphere 2018, 2018, 9, 9, 116 x FOR PEER REVIEW

18 of of 25 25 18

Figure 17. The stability parameter ζ distribution for the B01 (a–c), B09 (d–f), and B13 (g–i) bora episodes Figure 17. The stability parameter ζ distribution for the B01 (a–c), B09 (d–f), and B13 (g–i) bora at three different vertical levels: 2 m (a,d,g), 5 m (b,e,h), and 10 m (c,f,i). episodes at three different vertical levels: 2 m (a,d,g), 5 m (b,e,h), and 10 m (c,f,i).

The largest number of stability parameter values is grouped around zero, in accordance with The largest number of stability parameter values is grouped around zero, in accordance with the near-neutral thermal stratification of bora due to intensive mechanical mixing. This was also the near-neutral thermal stratification of bora due to intensive mechanical mixing. This was also found for a summer bora episode in [35] and is especially visible at two lower levels, where most found for a summer bora episode in [35] and is especially visible at two lower levels, where most of of the ζ are between −0.02 and 0.02. At the 10 m level, the frequency of those quasi-neutral cases the ζ are between −0.02 and 0.02. At the 10 m level, the frequency of those quasi-neutral cases is is lower compared to the 2 and 5 m level distributions, because more statically stable cases appear. lower compared to the 2 and 5 m level distributions, because more statically stable cases appear. Stable cases occur during nights when the heat flux is negative and the dynamical effects are either not Stable cases occur during nights when the heat flux is negative and the dynamical effects are either well developed or weakened [54]. The mean values of stability parameters are in general larger than not well developed or weakened [54]. The mean values of stability parameters are in general larger the corresponding medians for the selected bora episodes (Table 4). No significant difference related to than the corresponding medians for the selected bora episodes (Table 4). No significant difference the bora type was observed in the stability parameter distribution. related to the bora type was observed in the stability parameter distribution. Table 4. The summary statistics of the stability parameter ζ for each bora episode at selected Table 4.levels. The summary statistics of the stability parameter ζ for each bora episode at selected vertical vertical levels. Episode Episode HeightHeight (m) MeanMean Med (m) Med 2

B01 (SC) 5

B01 (SC)

10

2 5 B09 (DA) B09 (DA) 10 2 5 B13 (SA) B13 (SA) 10

2 5 10 2 5 10 2 5 10

0.0030.003 0.003 0.003 0.008 0.008 0.008 0.008 0.025 0.020

0.025

0.020

0.198

0.045

0.099

0.028

0.007 0.003 0.007 0.003 0.018 0.001 0.001 0.1980.018 0.045 0.027 0.004 0.004 0.0390.027 0.009 0.0990.039 0.028 0.009

Max Max 0.011 0.011 0.024 0.024 0.147 0.147 0.026 0.026 0.126 0.126 1.160 1.160 0.617 0.617 0.401 0.815 0.401 0.815

Min MinSD 0.0000.000 0.002 0.0000.000 0.005 −0.010 −0.010 0.028 0.000 0.000 0.008 −0.003 −0.003 0.038 −0.004 −0.004 0.307 −0.021 −0.021 0.106 −0.017 −0.301 −0.017 0.093 −0.301 0.223

SD 0.002 0.005 0.028 0.008 0.038 0.307 0.106 0.093 0.223

3.4. Monin–Obukhov Similarity Functions 3.4. Monin–Obukhov Similarity Functions The experimental values of the similarity functions (Φm, exp ) and the values obtained using the The experimental values of the similarity (𝛷in ) and18 thewith values obtained the 𝑚,𝑒𝑥𝑝 theoretical Equations (5) and (6) in Section 2.7functions are shown Figure respect to theusing stability theoretical Equations (5) and (6) in Section 2.7 are shown in Figure 18 with respect to the stability parameter ζ. parameter 𝜁.

Atmosphere 2018, 9, 116 Atmosphere 2018, 9, x FOR PEER REVIEW

19 of 25 19 of 25

Figure 18. The experimental similarity function function values values averaged The blue, blue, Figure 18. The experimental similarity averaged over over 30-min 30-min intervals. intervals. The magenta, and black markers denote the intervals from episode B01, B09, and B13, respectively, magenta, and black markers denote the intervals from episode B01, B09, and B13, respectively, compared to similarity functions (red(red lines), for the unstable (a,c) compared to the thetheoretical theoreticalvalues valuesofofthe the similarity functions lines), for statically the statically unstable and stable (b,d) surface layer. (a,c) and stable (b,d) surface layer.

The intervals stratification areare shown in the on the The intervalswith withunstable unstablethermal thermal stratification shown in diagrams the diagrams onleft, theand left,those and with stable stratification on the right. There is a significant discrepancy between the experimental those with stable stratification on the right. There is a significant discrepancy between and the the theoretical values, especially for ζ > 0 in the higher layer (5–10 m layer; Figure 18d). Furthermore, 2 experimental and the theoretical values, especially for 𝜁2 > 0 in the higher layer (5–10 m layer; the diagrams on the right (Figure 18b,d)on showing stable intervals relatively dispersion Figure 18d). Furthermore, the diagrams the right (Figure 18b,d)show showing stablehigh intervals show of the experimental values, especially in the higher level. Unfortunately, the diagrams showing the relatively high dispersion of the experimental values, especially in the higher level. Unfortunately, unstable intervals do notthe contain enough data points dispersion evaluation. Nevertheless, it can the diagrams showing unstable intervals do notforcontain enough data points for dispersion be noticed that the experimental values of the similarity functions are generally close to 1. This is to evaluation. Nevertheless, it can be noticed that the experimental values of the similarity functions be expected, since although the stability parameter is positive or negative, its absolute value stays are generally close to 1. This is to be expected, since although the stability parameter is positive or low—meaning that the stratification is mostly close neutral. negative, its absolute value stays low—meaning thattothe stratification is mostly close to neutral. In the lower layer (2–5 m), the experimental values have lower lower dispersion dispersion (Figure (Figure 18a,b). 18a,b). In the lower layer (2–5 m), the experimental values have Furthermore, they show a reasonable trend of increase when the lowest part of the atmosphere Furthermore, they show a reasonable trend of increase when the lowest part of the atmosphere becomes decrease when becomes more more unstable, unstable, which is in in accordance accordance becomes more more statically statically stable, stable, and and decrease when it it becomes which is with the theoretical expressions. This trend is not so obvious in the upper layer (5–10 m). It can can be with the theoretical expressions. This trend is not so obvious in the upper layer (5–10 m). It be concluded that the universal similarity functions are relatively successful at describing the wind profile concluded that the universal similarity functions are relatively successful at describing the wind in the lower of the butlayer, fail tobut give reliable above a certain Therefore, profile in thepart lower partsurface of the layer, surface fail to giveresults reliable results aboveheight. a certain height. the use of the universal similarity functions should probably be avoided in case of bora wind at Therefore, the use of the universal similarity functions should probably be avoided in case of or bora least windused or atwith least caution. used with caution. The reason functions do not givegive goodgood results in the in case of case bora wind most probably The reasonthe thesimilarity similarity functions do not results the of bora wind most lies in the failure of the main assumptions on which the similarity theory is based—that the surface probably lies in the failure of the main assumptions on which the similarity theory is based—that the layer is quasi-stationary, horizontal, and homogeneous. For a more would be surface layer is quasi-stationary, horizontal, and homogeneous. Fordetailed a moreanalysis, detailed itanalysis, it necessary to look at a larger number of bora episodes. would be necessary to look at a larger number of bora episodes. 3.5. Turbulence Kinetic Energy Budget 3.5. Turbulence Kinetic Energy Budget Figure 19 represents the TKE budget terms as defined in Equation (8) (Section 2.8) for the three Figure 19 represents the TKE budget terms as defined in Equation (8) (Section 2.8) for the three selected bora episodes: B01 (SC; Figure 19a), B09 (DA; Figure 19b) and B13 (SA; Figure 19c). selected bora episodes: B01 (SC; Figure 19a), B09 (DA; Figure 19b) and B13 (SA; Figure 19c).

Atmosphere 2018, 9, 116 Atmosphere 2018, 9, x FOR PEER REVIEW

20 of 25 20 of 25

(a)

(b)

(c) Figure 19. The turbulence kinetic energy (TKE) terms calculated on the two middle levels: 3.5 m (left) Figure 19. The turbulence kinetic energy (TKE) terms calculated on the two middle levels: 3.5 m (left) and 7.5 m (right). Term III (mechanical production) in the solid black line, term II (buoyant and 7.5 m (right). Term III (mechanical production) in the solid black line, term II (buoyant production) production) is the solid blue, dissipation is the solid red line, and the residual term is in dashed red. is the solid blue, dissipation is the solid red line, and the residual term is in dashed red. (a) B01 (SC), (a) B09 B01 (DA), (SC), and (b) B09 (DA), andThe (c)correlation B13 (SA). The correlation coefficient between production the mechanical (b) (c) B13 (SA). coefficient between the mechanical and production and dissipation is ~−0.9, and between the buoyancy term and residual is 0.8 3.5 andand 0.5 7.5 on dissipation is ~−0.9, and between the buoyancy term and residual is 0.8 and 0.5 on the the 3.5 and 7.5respectively. middle levels, respectively. middle levels, 2· 2 ·s − 3 during The shear shear term termdominates dominatesin inall allthree threeepisodes episodes(with (witha amaximum maximummagnitude magnitudeofof1 1mm s−3 during The the most intense phase of the selected bora episodes) meaning that the kinetic energy is extracted the most intense phase of the selected bora episodes) meaning that the kinetic energy is extracted from the the mean mean flow flow and and transformed transformed into into the the TKE. TKE. To To preserve preserve the the TKE TKE balance, balance, the the residual residual term term from 22 −3− 3 (with a minimum magnitude of −0.3 m · s ) and dissipation term (with a minimum magnitude of −1 (with a minimum magnitude of −0.3 m ·s ) and dissipation term (with a minimum magnitude of m2·s−3) are mostly negative for all episodes, decreasing the TKE in the layer considered here. Negative mechanical production (and the positive residual term), which is visible in the B09 episode,


21 of 25

−1 m2 ·s−3 ) are mostly negative for Atmosphere 2018, 9, x FOR PEER REVIEW

all episodes, decreasing the TKE in the layer considered here. 21 of 25 Negative mechanical production (and the positive residual term), which is visible in the B09 episode, is probably due to of u’ u’ in in the 30-min interval interval [55]. [55]. We We also is probably due to the the local local non-stationarity non-stationarity of the corresponding corresponding 30-min also found a minimum in the TKE (not shown) and u time series (Figure 14) corresponding to this block ∗ found a minimum in the TKE (not shown) and 𝑢∗ time series (Figure 14) corresponding to this block interval with interval with aa negative negative mechanical mechanical production. production. The other TKE terms together contribute with only only aa small (20%) to The other TKE terms together contribute with small portion portion (20%) to balancing balancing the the TKE TKE equation. Terms II (the buoyant production/consumption) and IV (vertical turbulent transport) are equation. Terms II (the buoyant production/consumption) and IV (vertical turbulent transport) are a a feworders ordersofofmagnitude magnitudesmaller smallerthan thanother otherTKE TKEbudget budgetterms. terms. Term Term IIII is is dominantly dominantly negative, negative, few acting as acting as aa weak weak sink, sink, while while term term IV IV is is dominantly dominantly positive. positive. These results agree well with the results from [49][49] for afor Mountain-Wave event event duringduring the T-REX These results agree well with the results from a Mountain-Wave the experiment. A newAoutcome is that is the terms only vary magnitude depending on theon bora T-REX experiment. new outcome that the terms onlyin vary in magnitude depending the type, bora while the signs of theof terms do notdo depend on theon synoptic type. type. type, while the signs the terms not depend the synoptic Comparing the two middle levels, absolute values of TKE terms on the 3.5 m, level Comparing the two middle levels, absolute values of TKE terms on lower, the lower, 3.5middle m, middle are larger than the corresponding 7.5 m middle-level values. This is due facttothat level are larger than the corresponding 7.5 m middle-level values. Thistoisthe due the turbulent fact that motions are more intense near ground level. The most intense turbulent motions are for B01—shallow turbulent motions are more intense near ground level. The most intense turbulent motions are for cyclonic bora type—followed by B09—deep by anticyclonic bora case. B01—shallow cyclonic bora type—followed B09—deep anticyclonic bora case. It is interesting to inspect the temporal correlation coefficients between mechanical production It is interesting to inspect the temporal correlation coefficients the between the mechanical term and dissipation, and between the buoyancy term and residual. The first correlation coefficient is production term and dissipation, and between the buoyancy term and residual. The first correlation large and negative (~−0.9 on both middle levels), whichlevels), can also be seen Figure 19. Furthermore, coefficient is large and negative (~−0.9 on both middle which canin also be seen in Figure 19. these two terms dominantly balance the TKE budget equation. This implies that the mechanical Furthermore, these two terms dominantly balance the TKE budget equation. This implies that the production dissipation aredissipation of a similar are size of anda are reciprocal. The second correlation coefficient mechanical and production and similar size and are reciprocal. The second mentioned is positive (0.8 and 0.5 on the 3.5 m and 7.5 m middle levels, respectively); consequently, correlation coefficient mentioned is positive (0.8 and 0.5 on the 3.5 m and 7.5 m middle levels, for statically stable conditions,for thestatically pressurestable transport is negative, for statically unstable conditions, respectively); consequently, conditions, the and pressure transport is negative, and it isstatically positive. unstable conditions, it is positive. for 3.6. Wind Profiles 3.6. Wind Profiles The vertical wind speed profiles reconstructed with the mean and median values (Section 2.9) are The vertical wind speed profiles reconstructed with the mean and median values (Section 2.9) shown in Figure 20. Both methods of reconstruction gave results with high (between 0.9874 are shown in Figure 20. Both methods of reconstruction gave results withcorrelations high correlations (between and 0.9995) and small relative errors. 0.9874 and 0.9995) and small relative errors.

Figure 20. The measured (circle) and reconstructed (line) vertical wind speed profiles in the Figure 20. The measured (circle) and reconstructed (line) vertical wind speed profiles in the x-direction x-direction with a percentage of near-neutral 30-min intervals, relative errors, and correlations. The with a percentage of near-neutral 30-min intervals, relative errors, and correlations. The profiles profiles reconstructed with the mean (median) values of 𝑢∗𝑝 (friction velocity estimated from wind reconstructed with the mean (median) values of u∗ p (friction velocity estimated from wind profile) and profile) and 𝑧0 are blue (green). The relative errors and correlations between the measurements and z0 are blue (green). The relative errors and correlations between the measurements and reconstructions reconstructions are given in the corresponding colors. are given in the corresponding colors.

In the examples depicted, the reconstructions of the vertical wind speed profiles based on the mean values have slightly better results, but that is not the general case. The efficiency of the method depends on the particular bora episode considered. This seems to be a consequence of wind speed,


22 of 25

In the examples depicted, the reconstructions of the vertical wind speed profiles based on the mean values have slightly better results, but that is not the general case. The efficiency of the method depends on the particular bora episode considered. This seems to be a consequence of wind speed, rather than bora type (different flow depth dynamics). When the vertical profiles shown in Figure 20 are normalized by the maximum average wind speed (not shown), they have the same shape, regardless of the bora type. This does not confirm the findings from [35], where different vertical profile shapes were observed for different wind speeds. We also tried reconstructing the wind profiles with the measured values of u∗ (calculated from turbulent fluxes using Equation (3)), but since these values are persistently lower than the values of u∗ p , such profiles underestimate the wind speed at all levels (not shown). 4. Conclusions We carried out, for the first time, a detailed analysis of high-frequency (20 Hz, downsampled to 10 Hz) wind data for several bora episodes measured at the Maslenica Bridge site in Croatia during autumn and winter 2015/2016 on three vertical levels (2, 5, and 10 m). A total of 14 bora episodes were detected and classified by depth and synoptic type, of which three typical episodes were selected and presented in this study: B01 (shallow cyclonic), B09 (deep anticyclonic), and B13 (shallow anticyclonic). Our results confirm the majority of the previous results [28,33,35]. The minimum energy (energy “gap”) in the frequency-weighted power spectral density graphs for the majority of the episodes is located at 30 min. We could not find clear evidence that it depends only on the bora type. Furthermore, power spectral densities disclose energy peaks at periods between 2 and 8 min for all three episodes, which are most likely related to bora pulsations. This further implies that mountain wave breaking occurred in all analyzed episodes. The thermal stratification during a bora episode is near-neutral due to intensive mechanical mixing, independent of the type of the episode. Deviations from this can be seen at the 10 m height in the nighttime, when the most statically stable bora cases occur. However, these never go beyond weak stratifications. The use of similarity functions in the bora surface layer was also tested. We suggest adopting the similarity theory for bora episodes with caution, since they fail to give reliable results, especially above a certain height. This is probably due to the fact that the main assumptions of the similarity theory are violated (i.e., quasi-stationarity and horizontal homogeneity). The vertical wind speed profiles—reconstructed with mean and median values—agree well with the logarithmic profile for the surface layer during all analyzed bora episodes. In the TKE equation, the shear term dominates in all three episodes extracting the kinetic energy from the mean flow and transforming into the TKE. The shear term is mainly balanced by the pressure transport (residual) and dissipation term. In the small set of typical bora types we analyzed (SC, DA, and SA), we found no evidence that possible differences in micro-scale properties are related to different bora types—which was one of the main goals of the study. The inspected elements that explicitly depend on the wind speed (i.e., friction velocity, TKE, and vertical wind profiles) are different, but that is not necessarily a function of bora type. The friction velocity and TKE budget terms increase with the increase of the mean streamwise wind component. In this paper, we present a novel approach to bora time series analysis, but we are aware of the possible limitations in finding micro-scale differences for different flow depth dynamics. The very sparse time series of sounding data (00 UTC and 12 UTC only) and the one-point measurements at Maslenica are a few of them. In order to further enhance this study (e.g., in the application of the similarity function), future work should aim for a more precise flow depth analysis, and perhaps more complex classification (e.g., by considering vertical wind shear and stability); but above all, more cases and measurements in a denser grid are needed to account for horizontal inhomogeneity. Furthermore, this study showed that the strongest bora episodes are mainly transitional in type (possible change in


23 of 25

flow depth dynamics), so to investigate the micro-scale properties of such cases, episodes should be divided into parts according to flow depth and synoptic type, and then analyzed. Acknowledgments: V.Š. and J.S. acknowledge Croatia Control Ltd. for their support. The work of A.B. has been supported by the Croatian Science Foundation (VITCLIC) project (PKP-2016-06-2975) which is funded by the Environmental Protection and Energy Efficiency Fund under the Government Program (Ministry of Environment and Energy & Ministry of Science and Education) for the Promotion of Research and Development Activities in the Field of Climate Change for the period 2015–2016. K.Š. acknowledges the EKONERG for their support. I.N., S.B., J.S., and M.B. acknowledge the Meteorological and Hydrological Service for their support. This work started as a project assignment within turbulence course at the Department of Geophysics, Faculty of Science, University of Zagreb. The study was partly funded by the Croatian Science Foundation projects CATURBO, No. 09/151 and HRZZ-IP-2016-06-2017 (WESLO). Author Contributions: All authors made a significant contribution to the data analysis, writing, and revising the manuscript. Vinko Šoljan coordinated the work, while also doing data analysis and writing. Andreina Belušić was responsible for the TKE analysis and writing. Kristina Šarović tested the Monin–Obukhov similarity theory. Irena Nimac analyzed the surface layer stability. Stjepana Brzaj analyzed the wind profiles. Jurica Suhin was responsible for the spectral analysis. Martin Belavić performed the data cleaning and quality control. Željko Veˇcenaj and Branko Grisogono designed and led the project, performed the wind data collection, and contributed to writing the manuscript. Conflicts of Interest: The authors declare no conflicts of interest.

References 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14.

15. 16. 17. 18.

Lepri, P.; Vecenaj, Z.; Kozmar, H.; Grisogono, B. Bora wind characteristics for engineering applications. Wind Struct. 2017, 24, 579–611. [CrossRef] Smith, R.B. Aerial observations of the Yugoslavian bora. J. Atmos. Sci. 1987, 44, 269–297. [CrossRef] Defant, F. Local winds. Compendium of Meteorology; American Meteorological Society: Boston, MA, USA, 1951; pp. 655–672. Jurˇcec, V. On mesoscale characteristics of bora conditions in Yugoslavia. In Weather and Weather Maps; Birkhäuser: Basel, Switzerland, 1981; pp. 640–657. Kuettner, J.; O’ Neill, T. ALPEX-the GARP mountain subprogram. Bull. Am. Meteorol. Soc. 1981, 62, 793–805. Smith, R.B. On severe downslope winds. J. Atmos. Sci. 1985, 42, 2597–2603. [CrossRef] Durran, D.R. Another look at downslope windstorms. Part I: The development of analogs to supercritical flow in an infinitely deep, continuously stratified fluid. J. Atmos. Sci. 1986, 43, 2527–2543. [CrossRef] ´ Bencetić Klaić, Z.; Belušić, D.; Grubišić, V.; Gabela, L.; Coso, L. Mesoscale airflow structure over the northern Croatian coast during MAP IOP—A major bora event. Geofizika 2003, 20, 23–61. Grubišíc, V. Bora-driven potential vorticity banners over the Adriatic. Q. J. R. Meteorol. Soc. 2004, 130, 2571–2603. [CrossRef] Gohm, A.; Mayr, G. Numerical and observational case-study of a deep Adriatic Bora. Q. J. R. Meteorol. Soc. 2005, 131, 1363–1392. [CrossRef] Grisogono, B.; Belušić, D. A review of recent advances in understanding the meso-and microscale properties of the severe Bora wind. Tellus A 2009, 61, 1–16. [CrossRef] Jurˇcec, V. The Adriatic frontal bora type. Croat. Meteorol. J. 1988, 23, 13–25. Kuzmić, M.; Grisogono, B.; Li, X.; Lehner, S. Examining deep and shallow Adriatic bora events. Q. J. R. Meteorol. Soc. 2015, 141, 3434–3438. [CrossRef] Lukšić, I. Tipovi strujanja zraka iznad Zagreba za vrijeme bure na Sjevernom Jadranu (Upper Level Flow Types at Zagreb During the Bora in the Northern Adriatic). In Papers from VII Climatological Conference; Federal Hydrometorological Institute: Belgrade, Serbia, 1972; pp. 111–129. Gohm, A.; Mayr, G.J.; Fix, A.; Giez, A. On the onset of bora and the formation of rotors and jumps near a mountain gap. Q. J. R. Meteorol. Soc. 2008, 134, 21–46. [CrossRef] Pullen, J.; Doyle, J.D.; Haack, T.; Dorman, C.; Signell, R.P.; Lee, C.M. Bora event variability and the role of air-sea feedback. J. Geophys. Res. Oceans 2007, 112. [CrossRef] Poje, D. Wind persistence in Croatia. Int. J. Climatol. 1992, 12, 569–586. [CrossRef] Pasarić, M.; Orlić, M. Djelovanje atmosfere na Jadran: danas te u uvjetima oˇcekivanih klimatskih promjena. Geofizika 2004, 21, 69–87.


19.

20. 21. 22. 23. 24. 25.

26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38.

39. 40. 41.

42.

43.

24 of 25

Belušić, A.; Prtenjak, M.T.; Güttler, I.; Ban, N.; Leutwyler, D.; Schär, C. Near-surface wind variability over the broader Adriatic region: Insights from an ensemble of regional climate models. Clim. Dyn. 2017, 1–26. [CrossRef] Mohorovicic, A. Interesting cloud pictures over the Bay of Buccari (with a comment from the editor J. Hann) (Interessante Wolkenbildung ber der Bucht von Buccari). Meteorol. Z 1889, 24, 56–58. Yoshino, M. Local Wind Bora; University of Tokyo Press: Tokyo, Japan, 1976. Stiperski, I.; Ivanˇcan-Picek, B.; Grubišić, V.; Bajić, A. Complex bora flow in the lee of Southern Velebit. Q. J. R. Meteorol. Soc. 2012, 138, 1490–1506. [CrossRef] Glasnović, D.; Jurˇcec, V. Determination of upstream bora layer depth. Meteorol. Atmos. Phys. 1990, 43, 137–144. [CrossRef] Klemp, J.B.; Durran, D.R. Numerical modelling of Bora winds. Meteorol. Atmos. Phys. 1987, 36, 215–227. [CrossRef] Watanabe, K. Bora and Man: Weather Forecasts and Prognosis of Bora by the Fishermen’s Traditional Way of Observation on the Croatian Coast. In Local Wind Bora; Yoshino, M.M., Ed.; University of Tokyo Press: Tokyo, Japan, 1976; pp. 267–274. Petkovšek, Z. Gravity waves and bora gusts. Ann. Meteorol. 1982, 19, 108–110. Petkovšek, Z. Main bora gusts—A model explanation. Geofizika 1987, 4, 41–50. Belušić, D.; Pasarić, M.; Orlić, M. Quasi-periodic bora gusts related to the structure of the troposphere. Q. J. R. Meteorol. Soc. 2004, 130, 1103–1121. [CrossRef] Belušić, D.; Žagar, M.; Grisogono, B. Numerical simulation of pulsations in the bora wind. Q. J. R. Meteorol. Soc. 2007, 133, 1371–1388. [CrossRef] Babić, N.; Veˇcenaj, Ž.; Kozmar, H.; Horvath, K.; De Wekker, S F.; Grisogono, B. On turbulent fluxes during strong winter bora wind events. Bound. Layer Meteorol. 2016, 158, 331–350. [CrossRef] Peltier, W.R.; Scinocca, J.F. The origin of severe downslope windstorm pulsations. J. Atmos. Sci. 1990, 47, 2853–2870. [CrossRef] Belušić, D.; Grisogono, B. Disappearance of pulsations in severe downslope windstorms. Croat. Meteorol. J. 2005, 40, 84–87. Veˇcenaj, Ž.; Belušić, D.; Grisogono, B. Characteristics of the near-surface turbulence during a bora event. Annales Geophysicae 2010, 28, 155–163. [CrossRef] Veˇcenaj, Ž.; Belušić, D.; Grubišić, V.; Grisogono, B. Along-coast features of bora-related turbulence. Bound. Layer Meteorol. 2012, 143, 527–545. [CrossRef] Lepri, P.; Kozmar, H.; Veˇcenaj, Ž.; Grisogono, B. A summertime near-ground velocity profile of the Bora wind. Wind Struct. 2014, 19, 505–522. [CrossRef] Lepri, P.; Veˇcenaj, Ž.; Kozmar, H.; Grisogono, B. Near-ground turbulence of the Bora wind in summertime. J. Wind Eng. Ind. Aerodyn. 2015, 147, 345–357. [CrossRef] Vickers, D.; Mahrt, L. Quality control and flux sampling problems for tower and aircraft data. J. Atmos. Ocean. Technol. 1997, 14, 512–526. [CrossRef] Ivanˇcan-Picek, B.; Grubišić, V.; Stiperski, I.; Xiao, M.; Bajić, A. “Zadar calm” during severe Bora. In Proceedings of the 29th International Conference on Alpine Meteorology, Extended Abstracts, Chambery, France, 4–8 June 2007. Telišman Prtenjak, M.; Belušić, D. Formation of reversed lee flow over the north-eastern Adriatic during bora. Geofizika 2009, 26, 145–155. Deutscher Wetterdienst (DWD) Surface Analysis Charts. Available online: http://www1.wetter3.de/ Archiv/archiv_dwd.html (accessed on 28 November 2016). National Centers for Environmental Prediction, National Weather Service and U.S. Department of Commerce. NOAA NCEP GDAS/FNL 0.25 Degree Global Tropospheric Analyses and Forecast Grids. Available online: https://doi.org/10.5065/D65Q4T4Z (accessed on 10 March 2017). University of Wyoming Department of Atmospheric Science Sounding Data from University of Wyoming Database. Available online: http://weather.uwyo.edu/upperair/sounding.html (accessed on 28 November 2016). Reynolds, O. On the dynamical theory of incompressible viscous fluids and the determination of the criterion. Proc. R. Soc. Lond. 1894, 56, 40–45. [CrossRef]


44.

45. 46. 47.

48. 49. 50. 51. 52. 53. 54. 55.

25 of 25

Stull, R.B. An Introduction to Boundary Layer Meteorology, 1st ed.; Stull, R.B., Ed.; Atmospheric and Oceanographic Sciences Library; Springer: Dordrecht, The Netherlands, 1988; Volume 13, ISBN 978-90-277-2768-8. Welch, P. The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms. IEEE Trans. Audio Electroacoust. 1967, 15, 70–73. [CrossRef] Kaimal, J.C.; Finnigan, J.J. Atmospheric Boundary Layer Flows: Their Structure and Measurement; Oxford University Press: New York, NY, USA, 1994. Albertson, J.D.; Parlange, M.B.; Kiely, G.; Eichinger, W.E. The average dissipation rate of turbulent kinetic energy in the neutral and unstable atmospheric surface layer. J. Geophys. Res. Atmos. 1997, 102, 13423–13432. [CrossRef] Fortuniak, K.; Pawlak, W.; Siedlecki, M. Integral turbulence statistics over a central European city centre. Bound. Layer Meteorol. 2013, 146, 257–276. [CrossRef] Veˇcenaj, Ž.; De Wekker, S.F.; Grubišić, V. Near-surface characteristics of the turbulence structure during a mountain-wave event. J. Appl. Meteorol. Climatol. 2011, 50, 1088–1106. [CrossRef] Tennekes, H.; Lumley, J.L. A First Course in Turbulence; MIT Press: Cambridge, MA, USA, 1972. Holton, J.R. An Introduction to Dynamic Meteorology, 4th ed.; Dmowska, R., Holton, J.R., Eds.; International Geophysics Series; Elsevier Academic Press: Burlington, MA, USA, 2004; ISBN 978-0-12-354015-7. Babić, K.; Bencetić Klaić, Z.; Veˇcenaj, Ž. Determining a turbulence averaging time scale by Fourier analysis for the nocturnal boundary layer. Geofizika 2012, 29, 35–51. Belušić, D.; Pasarić, M.; Pasarić, Z.; Orlić, M.; Grisogono, B. A note on local and non-local properties of turbulence in the bora flow. Meteorol. Z. 2006, 15, 301–306. [CrossRef] Tampieri, F. Turbulence and Dispersion in the Planetary Boundary Layer; Physics of Earth and Space Environments; Springer International Publishing: Cham, Switzerland, 2017; ISBN 978-3-319-43602-9. Veˇcenaj, Ž.; De Wekker, S.F. Determination of non-stationarity in the surface layer during the T-REX experiment. Q. J. R. Meteorol. Soc. 2015, 141, 1560–1571. [CrossRef] © 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).