Investigation of an extreme Koshava wind episode of 30 January–4 February 2014

An extreme Koshava episode (EKE) from 30 January to 4 February 2014 has been studied. Koshava is a local windstorm in Southeast Europe. EKE was characterized by wind gusts above 45 m s−1 and deep snowdrifts. Strong Eurasian anticyclone in combination with large temperature gradient between the anticyclone region and the Mediterranean caused the EKE. The anticyclone had the probability of occurrence of 0.1%. Koshava layer was either statically stable or adiabatic. The event was numerically modelled using two mesoscale models: WRF‐NMM and NMMB. Wind directions were forecasted more accurately than the mean speed and gusts.


Introduction
From 30 January to 4 February 2014, large parts of Serbia were affected by an extremely severe windstorm. A strong wind, named Koshava, occasionally had hurricane velocities. In addition to extreme wind speeds, Koshava created snowdrifts several meters deep. This weather disaster impacted everyday life across the country -trees were toppled, many buildings and cars were damaged, and rail and air transports in northern Serbia were shut down. Furthermore, between 5000 and 10 000 people were trapped in their vehicles on the snow-covered roads and it was necessary to rescue them from cold. Fortunately, no human casualties were reported. Interestingly, this extreme Koshava episode (EKE) had a positive effect on the air quality in cities. Measurements in Belgrade showed that the atosphere after the EKE was two times cleaner compared with the air quality before the event.
Koshava is a strong windstorm that blows from southeast directions over most regions of Serbia. The southeasterly winds need to last for at least 2 days and possess the mean hourly wind speeds above 2 m s −1 to be considered as Koshava (Romanić et al., 2015a). Romanić et al. (2015b) demonstrated that the Eurasian anticyclones and Mediterranean cyclones, together with the orography of the Balkan region are the main Koshava drivers. The pressure gradients created by these two pressure systems move air masses from Ukraine and Moldavia towards the Adriatic Sea. However, the South Carpathian and Balkan Mountains redirect that air towards the Iron Gates ( Figure 1). The diverging jets downstream from the Iron Gates are known as the Koshava wind.
Koshava is strongest in the Vršac (VR) region. Figure 2 shows that the VR station has 10-15 m s −1 stronger Koshava gusts than the Belgrade (BG) station. The strongest Koshava gust ever recorded was measured in VR on 11 February 1987 [48 m s −1 at 10 m above ground (a.g)]. The strongest Koshava speed at the BG station was 35.9 m s −1 at 24 m a.g. on 17 October 1976. The measured velocity at 650 m a.g. was 45 m s −1 . These two weather stations mostly recorded different storms as being the most extreme. Exceptions are the storms that occurred in February 1979 and February 2014. These results indicate that EKEs are mostly localized. An absence of extreme Koshava winds in the period from the late 1980s to 2014 is in accordance with the results by Romanić et al. (2015a), which state that Koshava was most active in the period from the 1970s to 1980s.
This article investigates meteorological factors that caused the EKE in late January to early February of 2014. This EKE is the most extreme wind event in Serbia that has been peer-reviewed analysed. The following questions are addressed: (1) what was the synoptic situation which brought EKE, (2) how well two mesoscale models with different configurations forecasted the EKE, and most importantly, (3) what caused such high Koshava speeds to occur? Hereafter, when refereeing to dates related to the EKE, the year (2014) will be omitted for simplicity.   1970 1975 1980 1985 1990 1995 2000 2005 2010 2015 Wind speed (m s -1 ) Year Figure 2. Past records of the extreme Koshava gusts at Vršac station (blue diamonds) and Belgrade station (red stars) Serbia (RHMSS). The mean sea level pressure (MSLP) and 2-m air temperature fields were obtained from the European Center for Medium-Range Weather Forecasting's (ECMWF) Meteorological Archival and Retrieval System database.
This study uses the Weather Research and Forecasting (WRF) model for the numerical simulation of the EKE. The employed solver is the Nonhydrostatic Mesoscale Model (NMM; version 3.6.1), which solves the compressible, nonhydrostatic Navier-Stokes equations on Arakawa-E grid. The vertical coordinate is a hybrid sigma-pressure coordinate. Two different configurations of the NMM model were tested (Table 1). In addition, two more numerical simulations were performed using a modified version of the NMM model on Arakawa-B grid, called NMMB. In total, four simulations with different models' configurations are performed -two using NMM and two using NMMB. The forecasts (mean hourly wind speed, wind gust, and wind direction) were compared against measurements from four weather stations situated in the Koshava region (BG, Novi Sad (NS), VR, and Veliko Gradište (VG)). The results are also benchmarked against the ECMWF's Integrated Forecast System (IFS) analysis.
Models' configurations are summarized in Table 1. NMM, NMMB_1, and NMMB_2 setups are operationally used at the RHMSS. Due stably stratified planetary boundary layer (PBL; see Section 3.1), the NMM_QNSE configuration uses the Quasi-Normal Scale Elimination scheme (QNSE; Sukoriansky et al., 2005) for PBL parametrization. The NMM and NMMB simulations have similar physics, but the forecasts are based on different initial and boundary conditions.

Microphysics
Ferrier (Ferrier et al., 2002) Longwave radiation GFDL (Schwarzkopf and Fels, 1991) RRTMG (Iacono et al., 2008) Shortwave radiation GFDL (Lacis and Hansen, 1974) RRTMG (Iacono et al., 2008 Namely, the quality of forcing data can have a pronounced influence on the accuracy of forecasts (Talbot et al., 2012). Numerical modelling of surface and near-surface variables under stable atmospheric conditions is more challenging compared with unstable stratification (Shin and Hong, 2011). This difficulty partially arises because turbulent eddies in stable atmosphere are much smaller compared with unstable (convective) conditions (Beare et al., 2006). For that reason, the NMM_QNSE simulation relies on the QNSE PBL scheme. This scheme is particularly designed for stable stratifications. Namely, QNSE is a turbulent kinetic energy prediction scheme based on a quasi-Gaussian spectral closure model with the surface layer parameters derived from large eddy simulations of stably stratified atmosphere (Sukoriansky et al., 2005).

Synoptics and dynamics of the EKE
Animation S1 (Supporting Information) shows the ridge of a deep anticyclone stretching across the Balkan Peninsula on 27 January. On 30 January (Figure 3a), the pressure difference between the Eurasian anticyclone and Mediterranean cyclone was 55 mb, and the temperature difference between these two regions was approximately 50 ∘ C. The distance between the centres of these two pressure systems was about 2500 km. The pressure gradients above Serbia were extremely large -approximately 5 mb/100 km. This synoptic situation was favourable for the development of a strong Koshava wind.
EKE started in the afternoon of 30 January. The central pressure in the anticyclone that caused the EKE was 1055 mb, which according to Romanić et al. (2015b) has an exceedance probability of occurrence of 0.1%. The Mediterranean cyclone, situated in the Gulf of Genoa, was not particularly strong as its central pressure was 3 mb below the average central pressure of the cyclones which typically generate Koshava (Romanić et al., 2015b). Therefore, the strong anticyclone was the main trigger for the EKE. Animation S1 further shows that while the centre of the anticyclone was slowly moving eastward from the Koshava region during the EKE, its high-pressure ridge was constantly positioned above Balkan Peninsula and therefore maintained strong pressure and temperature gradients above the region. The disappearance of the anticyclone (Figure 3b) led to the cessation of the EKE.
The Synoptic Koshava Index (SKI; Romanić et al. (2015b)), defined as the normalized area-averaged MSLP difference between the anticyclone and cyclone regions is presented in Figure 4a. The SKI values above 25 are noticeable during EKE while the values above 35 can be observed for the days with the strongest Koshava wind. The high values of SKI demonstrate a strong pressure difference between the anticyclone and cyclone regions, which indubitably generated the strong Koshava velocities. Romanić et al.'s (2015b) results show that the SKI values above 35 have a probability of occurrence below 1% (less than 0.3% chance of the SKIs above 40). This case study suggests that the SKI is a reliable prognostic Koshava parameter.
The major Koshava contributors on the mesoscales are the across-mountain MSLP and potential temperature differences between west parts of the Wallachia Valley and the Koshava region (Romanić et al., 2015b). The across-mountain MSLP difference between Drobeta-Turnu Severin and VG reached 10 mb during the EKE, whereas the across-mountain potential temperature difference between these two stations was below −7 ∘ C. The probability of occurrence of such a large pressure difference is 0.01%, while the joint probability of these values is practically zero (0.0012%; Romanić et al., 2015b).
Due to these large pressure and temperature differences, northern and northeastern parts of Serbia experienced a blizzard that caused significant damage to infrastructure and posed a threat to human safety. A strong and cold Koshava wind with daily mean speeds of 6-12 m s −1 and occasional gusts between  12 and 25 m s −1 was blowing during the EKE. In the early morning of 1 February, Koshava gusts in BG reached 29 m s −1 . Radiosonde measurements in BG conducted on 2 February (00 UTC) recorded the maximum Koshava speed of 34 m s −1 at a height of 497 m a.g (Figure 4b). The strongest Koshava speed, however, was measured in VR on 1 February, when a gust at 10 m a.g. reached 47 m s −1 . This value is the second largest Koshava gust ever recorded. The 2014 EKE was the first severe Koshava event after more than 25 years without any similar storm ( Figure 2). Emagrams in Figure 5 indicate that the Koshava layer at the BG station was considerably deeper than at the VR station. The Koshava winds above BG occurred in a layer stretching from the surface up to 350 mb. The Koshava layer at the VR station was much shallower; reaching only up to 900 mb. Despite the differences in thickness, the maximum Koshava speeds occurred at similar heights (950-900 mb) at both stations. Figure 4 shows that the maximum Koshava speed at the BG station was measured in a layer between 389 (31 January) and 594 m a.g. (3 February). This observation is in accordance with the finding by Vukmirović (1985) that Koshava is strongest in the layer between 500 and 600 m a.g. The strongest Koshava speed in BG on 2 February occurred at the bottom of an elevated inversion, which was also previously noticed by Unkašević et al. (1999). Koshava layer was characterized either with a deep surface inversion (Figure 5b) or an elevated inversion layer above adiabatic and conditionally stable layers close to the surface (Figure 5a). The elevated inversion was due to the advection of the warm air in the front of the Mediterranean cyclone (Animation S1). Atmospheric stability indices in emagrams show the absence of the atmospheric instability during the EKE. However, both environmental helicity (EH) and storm relative environmental helicity (SREH) had very high values. The helicity represents a measure of the vertical transfer of energy by the shear of the horizontal wind vector. The SREH values at the VR station were around 370 m 2 s −2 , which is comparable with the figures favourable for the development of strong tornadoes. Figure 6 and Table 2 show the models' forecasted wind directions (D) more precisely than the mean hourly . Parameters in box as follows: K, K index in ∘ C; TT, total totals index in ∘ C; PW, precipitable water for the entire sounding in centimetres; Temp, temperature on ground in ∘ C; Dewp -dewpoint on ground in ∘ C; Thetae, equivalent potential temperature in K; LI, lifted index in ∘ C; CAPE, convective available potential energy in J kg −1 ; CIN, convective inhibition in J kg −1 ; EH, environmental helicity in m 2 s −2 ; SREH, storm relative environmental helicity in m 2 s −2 ; StrmDir, storm direction in degrees; StrmSpd, storm speed in m s −1 . For further explanation of the parameters see Doswell III and Schultz (2006). wind speeds (V) and wind gusts (V). The most accurate wind direction forecasts were at the VR station. A constant 20-25 ∘ offset between the forecasted (∼130 ∘ ) and the observed (∼110 ∘ ) wind directions at the VG and NS stations are evident. The forecasts at BG were the least accurate, probably due to the complex urban environment around the BG station. The higher accuracy of wind direction forecasts compared with wind speed forecasts is a result of the strong pressure gradients that drove Koshava unidirectionally. The discrepancies between modelled and observed wind directions are most likely caused by poor resolution of small-scale orography and lack of proper representation of urban environments in the models.

Numerical modelling of the EKE
The models were not fully capable of capturing high fluctuations of the observed wind speeds and gusts. The most precise V forecasts were at the BG and NS stations, i.e. the stations that are further downstream from the mountain regions of central and eastern Serbia. A general tendency of NMM_QNSE to predict higher V andV values than the other four models is noticeable. Small variations between NMM and NMMB simulations indicate that initial and boundary conditions, as well as the grid resolution, had little impact on the forecasts' accuracy. The largest discrepancies between modelled and observed values of V are found at the VR station, where the observed Vs were approximately two times higher than the modelled values in the first 90 h of simulation. The largest deviations between NMM_QNSE's outputs and forecasts of the other four models were in cities (the BG and NS stations). NMM produced the most accurate forecasts of V at the VG station. The fluctuations of mean wind speed were best modelled at the BG and VR stations (high CC values). The superiority of the QNSE scheme for numerical modelling of stable conditions, however, has not been confirmed in this study. It can be concluded that numerical modelling of stable PBLs is not at a satisfactory level and more research is required.
The wind gusts are not directly calculated by WRF, but computed in post-processing. The Unified Post Processor (UPP; version 2.2) is used to post-process NMM and NMMB outputs. UPP computes surface wind gusts (V) by mixing down momentum from the PBL height (z PBL ) level (Mankin, 2015): ) .
Here, V PBL is the wind speed at z PBL and V sfc is the surface wind speed. The assumption behind this method is that the gusts are caused by air parcels brought down from the top of the PBL by turbulent eddies. The reliability of calculated gusts, therefore, depends on the accuracy of calculated turbulent kinetic energy and other PBL features. Applicability of the method is particularly questionable for the stable PBL (hence small turbulent eddies) and extremely large helicities, both observed during the EKE. However, analysis of  Table 1: NMM (green lines with stars), NMMB_1 (blue lines with triangles), NMMB_2 (red lines with diamonds), NMM_QNSE (cyan lines with dots) and IFS (black lines with circles). See Table 2 for for the verification statistics.
different methods for calculations of wind gusts is beyond the scope of this study (c.f. Sheridan, 2011).
RMSDs of estimated gusts are as high as 12.9 m s −1 (NMM_QNSE's forecasts at the VR stations). The discrepancies between observed and simulated gusts are smallest at the BG station were NMM_QNSE captured gusts fairly well. In other cases, gusts were largely either over-predicted (NS and VG) or under-predicted (VR). The presented results demonstrate that more research is needed in the field of wind gust modelling. Table 2. Bias (i.e. absolute difference), mean absolute difference (MAD), root mean square difference (RMSD) and correlation coefficient (CC) between forecasts and observations during the EKE. The corresponding time series are portrayed in Figure 6.    Figure 7 is spatial distribution of the mean daily wind speeds in the Koshava region based on the NMM simulations. The strongest Koshava speeds (>17 m s −1 ) occurred in east Serbia (the VR region) and west Romania (the region around Caransebes). High wind speeds were also found in the Morava River basin (14-21 m s −1 ). The high intensity of the Koshava wind in the above-mentioned areas is due to the wind channelling effect created by the local orography (Romanić et al., 2015b). Strong winds were also blowing along the gorges in the South Carpathians. Koshava weakened moving downstream from these regions, but strong winds were still present in around BG and NS. Extremely high Koshava speeds were occurring from 31 January to 2 February, when the mean daily wind speeds above the whole Pannonian Plane were greater than 10 m s −1 . Figure 7, however, should be interpreted cautiously as the RMSDs in Table 2 demonstrate that NMM greatly under-predicted Vs in the VR region.

Conclusions
This article is a case study of an EKE from 30 January to 4 February, 2014. The mean daily Koshava speeds during the EKE were generally above 10 m s −1 with gusts exceeding 30 m s −1 . Similar Koshava events had not occurred during the past 25 years. Meteorological factors leading to this EKE were: (1) large pressure gradients across the Balkan Peninsula caused by an extremely deep Eurasian anticyclone with the central pressure above 1055 mb combined with (2) large temperature gradients between the anticyclone and cyclone regions -all of which resulted in (3) extremely large across-mountain pressure and potential temperature differences between the Wallachia Valley and the Koshava region. The extreme wind speeds occurred together with deep snowdrifts, thus causing a blizzard. The atmosphere in the Koshava layer was either stable with surface inversion or adiabatic with capping elevated inversion.
Four numerical simulations of the EKE were performed (two using WRF-NMM and two using NMMB) and benchmarked against the observations and the global IFS' analysis. The models captured Koshava's directions more accurately than its mean speeds and gusts. The strongest Koshava speeds occurred in the Vršac region. Higher accuracy of direction forecasts is due to the strong pressure gradients that drove the wind. The low accuracy of modelled speeds and gusts is probably due to the lack of planetary boundary layer schemes capable of properly replicating stable stratifications. More research is needed in the field of numerical modelling of stable boundary layers and wind gusts.