Numerical modeling and analysis of the effect of complex Greek topography on tornadogenesis

Tornadoes have been reported in Greece over recent decades in specific sub-geographical areas and have been associated with strong synoptic forcing. While it has been established that meteorological conditions over Greece are affected at various scales by the significant variability of topography, the Ionian Sea to the west and the Aegean Sea to the east, there is still uncertainty regarding topography’s importance on tornadic generation and development. The aim of this study is to investigate the role of topography in significant tornadogenesis events that were triggered under strong synoptic scale forcing over Greece. Three tornado events that occurred over the last years in Thebes (Boeotia, 17 November 2007), Vrastema (Chalkidiki, 12 February 2010) and Vlychos (Lefkada, 20 September 2011) were selected for numerical experiments. These events were associated with synoptic scale forcing, while their intensities were T4–T5 (on the TORRO scale), causing significant damage. The simulations were performed using the non-hydrostatic weather research and forecasting model (WRF), initialized by European Centre for Medium-Range Weather Forecasts (ECMWF) gridded analyses, with telescoping nested grids that allow for the representation of atmospheric circulations ranging from the synoptic scale down to the mesoscale. In the experiments, the topography of the inner grid was modified by: (a) 0 % (actual topography) and (b) −100 % (without topography), making an effort to determine whether the occurrence of tornadoes – mainly identified by various severe weather instability indices – could be indicated by modifying topography. The principal instability variables employed consisted of the bulk Richardson number (BRN) shear, the energy helicity index (EHI), the storm-relative environmental helicity (SRH), and the maximum convective available potential energy (MCAPE, for parcels with maximumθe). Additionally, a model verification was conducted for every sensitivity experiment accompanied by analysis of the absolute vorticity budget. Numerical simulations revealed that the complex topography constituted an important factor during the 17 November 2007 and 12 February 2010 events, based on EHI, SRH, BRN, and MCAPE analyses. Conversely, topography around the 20 September 2011 event was characterized as the least significant factor based on EHI, SRH, BRN, and MCAPE analyses.


Introduction
Tornadoes are violently rotating columns of air, associated with a swirling cloud of debris and a funnel shaped cloud extending downward from the base of the parent cumulonimbus cloud.They are associated with extremely strong winds, inside and around the tornado's funnel, causing extended damage and in many cases, loss of life.Tornadoes and waterspouts have always fascinated mankind.They were well known to the ancients, as virtually all classical philosophers attempted explaining them.Similarly, over the last 40 years researchers have been engaged in studies and attempts to understand the mechanisms that develop such phenomena.
Published by Copernicus Publications on behalf of the European Geosciences Union.
Tornadogenesis is a major forecast challenge in regions of complex terrain (e.g., Homar et al., 2003).It is no wonder that the complex Greek inland terrain, along with the Ionian Sea to the west and the Aegean Sea to the east, appears to be a vulnerable area for convective weather and tornado occurrence.Operational and research meteorologists often refer to and investigate the presence of diagnostic variable sets that can be used as forecast parameters for severe convective weather.Diagnostic variables have a long history in relation to forecasting severe convective weather (Schaefer, 1986;Johns and Doswell, 1992) as their values represent their capacity to summarize in a single value some kinematic or thermodynamics characteristics of the severe storm environment.Doswell and Schultz (2006) have discussed clearly these diagnostic variables and introduced a five scale classification scheme: (1) simple observed variables, (2) simple calculated variables, (3) derivatives or integrals (spatial or temporal) of simple observed or calculated variables, (4) combined variables, and (5) indices.Are these diagnostic variables sufficient to forecast tornadogenesis?What other factors make the environment conditions favorable to tornadoes?
In order to answer these questions we need to review the diagnostic variables that have been identified by prior research as "ingredients" for tornado occurrence.The concept of ingredients was introduced by Doswell (1987), Doswell et al. (1996), and Groenemeijer et al. (2011), and can be defined as an essential, but not sufficient, condition for a phenomenon to occur.For example, a storm requires two ingredients: (1) convective available potential energy (CAPE; Table 1) and (2) lift.By lift, we mean the sufficient initial vertical motion that results in convection, and this could be carried out either by thermal processes or by a synoptic driven force (fronts) or even by the effect of topography/orography. Miglietta and Regano (2008) confirm the importance of even a relatively shallow orography for the development of convective cells, using weather research and forecasting model (WRF) sensitivity experiments over a small area in Apulia (in southern Italy).Regarding CAPE, several researchers have noted that large values of CAPE are not mandatory for tornado development (Monteverdi et al., 2003;Hannesen et al., 1998;Groenemeijer and Van Delden, 2007).Besides CAPE, Rasmussen (2003) has shown that the storm relative helicity (SRH) in the lower atmosphere (Table 1) can be used to distinguish non-tornadic from weakly and significantly tornadic storms.On the other hand, for Craven and Brooks (2002) and Brooks (2009) the strong low level wind shear was identified as a good predictor of tornadoes.Energy helicity index (EHI; Table 1) was evaluated by Rasmussen (2003) on the basis of three classes of proximity soundings, relating convective weather and supercell storms with tornadic and non tornadic events.The bulk Richardson number (BRN, Table 1) is used to quantify the relationship between buoyant energy and vertical wind shear (Moncrieff and Green, 1972), taking into account the wind components of the difference between the density-weighted mean winds over the 6000 m and the 500 m a.g.l.As discussed in Droegemeier et al. (1993), the BRN is only a gross estimate of the effects of vertical wind shear on convective storms, since it does not measure the turning of the wind profile with height.However, Weisman and Klemp (1984) show, using cloud-scale model simulations, that the BRN can distinguish between supercell and multicell storms, with modeled supercells likely when 10 ≤ BRN ≤ 50 and multicell storms likely when BRN > 35.
It is important to note that there is no well-defined threshold value for BRN, since there is an overlap in the values used to specify storm type.
A model simulation with spatial resolution that reproduces the weather conditions at storm scale is a powerful tool for understanding the relevant physical processes involved and calculating the above-mentioned diagnostic variables.In our study, three tornado events were simulated using the nonhydrostatic WRF model, in order to determine whether the model is able to indicate the occurrence of tornadic activity by modifying the topography.The selected tornado events have been associated with synoptic scale forcing; their intensities were T4-T5 on the TORRO scale (Meaden, 1976) and caused significant damage.The TORRO scale describes peak winds, yet it was designed as a formal extension of the B-scale (Meaden, 1976).T4 and T5 are characterized as significant/strong tornadoes with a wind velocity ranging between 51-61 m s −1 (T4) and 61-71 m s −1 (T5).The analyzed tornado events have occurred over recent years in three Greek areas; namely in Thebes (central Greece, Boeotia, 17 November 2007), Vrastema (north Greece, Chalkidiki, 12 February 2010) and Vlychos (west Greece, Lefkada, 20 September 2011).
An analysis of severe weather variables was conducted for every simulation.The diagnostic variables comprise the BRN, EHI, SRH, and maximum convective available potential energy (MCAPE, for parcels with maximum θ e ).The low-level SRH, can vary dramatically over short distances, due to topographically channeled flow (e.g., Braun and Monteverdi, 1991;Bosart et al., 2006).The selection of these parameters is consistent with the approach adopted in previous studies (Droegemeier et al., 1993;Johns et al., 1993;Rasmussen and Blanchard, 1998;Brooks et al., 2003;Doswell and Evans, 2003;Thompson et al., 2003;Shafer et al., 2009;Matsangouras et al., 2011aMatsangouras et al., , 2012;;Sioutas et al., 2012).(in hPa) this value is to be taken from; α denotes the specific volume and lp denotes a value associated with a lifted parcel; LFC stands for a lifted parcel's level of free convection and EL stands for its equilibrium level.For the Bulk Richardson number, U 0 denotes the density-weighted speed of the mean vector wind in the layer 0-6 km, and U denotes the speed of the mean vector wind in the layer from the surface to 500 m -the quantity is sometimes referred to as the "BRN shear".For the storm-relative helicity, k × dV h dz is the horizontal vorticity, and vector k is the unit vector in the vertical, V h − c is the storm-relative environmental wind velocity.The shear vector consists of two components, the speed shear V h (a change in wind speed with height) and the directional shear or storm motion vector (a change in wind direction with height).

Index name Formula Reference
Convective potential available energy (CAPE) CAPE = EL LFC (αlp − α) dp Glickman (2000, p. 176) Weisman and Klemp (1982) Storm relative helicity (SRH) Droegemeier et al. (1993), Markowski et al. (1998) Energy helicity index (EHI) EHI = (CAPE) (SRH)   160.000Hart and Korotky (1991) Overall, the goal of this paper is to investigate the importance of complex terrain in tornado formation.The data sources and WRF setup are presented in Sect. 2. Section 3 analyzes the selected tornado events and puts forwards a brief weather synoptic analysis.In Sect. 4 we show the model simulation verification.In Sect. 5 we demonstrate our analysis of absolute vorticity budget and in Sect.6 we examine and discuss the results obtained.Finally, Sect.7 summarizes our findings and conclusions.

Data sources and model setup
Starting in 2007, the Laboratory of Climatology and Atmospheric Environment (LACAE, http://lacae.geol.uoa.gr) at the University of Athens has made a systematic effort to record tornadoes, waterspouts, and funnel clouds over Greece.In 2009, LACAE developed an open-ended, online tornado reporting database web system (http://tornado.geol.uoa.gr), contributing to the compilation of a climatology concerning these extreme weather events.A flow chart of the LACAE tornado reporting data stream system, accompanied by a plausibility check process, was presented by Matsangouras et al. (2014).The Laboratory of Climatology and Atmospheric Environment is in close collaboration with the European Severe Storm Laboratory (ESSL) and regularly submits reports to the European Severe Weather Database (ESWD).
Tornado event simulations were performed using the nonhydrostatic WRF model with the dynamical core of advanced research WRF (ARW).An analysis of BRN, EHI, SRH, and MCAPE (for parcels with maximum θ e ) severe weather variables was conducted for every simulation.
The WRF-ARW V3.2 non-hydrostatic numerical model (Skamarock et al., 2008;Wang et al., 2010) was used in order to simulate tornadic activity weather conditions on 17 November 2007, 12 February 2010, and 20 September 2011.Three one-way nested domains were utilized.The spatial resolution of the model was 12 km for D1 (381×301 grid points), 4 km for D2 (310 × 301 grid points) and 1.333 km for D3 (202×202 grid points).The inner domain employed a very fine grid spacing in order to provide a more realistic representation of topography and atmospheric flow in the areas of interest, which are characterized by complex terrain.Mass et al. (2002) showed that increasing model resolution generally leads to improved simulations in regions with complex terrain.Similarly, the simulations of Colle andMass (1998, 2000a, b) were improved as the grid spacing decreased to 1-1.33 km and many of the mesoscale features were more realistically simulated.Figure 1 illustrates all nested domains for every case study.All nests were integrated in non-hydrostatic mode.In all cases domain D1 covered Europe and northern Africa and D2 covered an area from the southern Balkans (northern boundary) to the north coast of Africa (southern boundary), and from southern Italy to eastern Turkey.Domain D3 was setup over each area of interest.In the sensitivity experiments, topography was modified only in D3 domain according to the methodology of Miao et al. (2003) and Koletsis et al. (2009Koletsis et al. ( , 2010)).Chiao et al. (2004) and Chen et al. (2010Chen et al. ( , 2011Chen et al. ( , 2013) followed a similar procedure in their numerical experiments on the role of orography, with no modifications at the topography of their external grid.Figure 1 illustrates the actual topography of Greece and the three nested domains.In the vertical, 39 sigma levels (up to 50 hPa) with increased resolution in the boundary layer were utilized by all nests.ERA-Interim reanalysis data (Dee et al., 2011) from the European Centre for Medium-Range Weather Forecasts (ECMWF) with a spatial resolution of 0.75 • × 0.75 • were used as initial and lateral boundary conditions for the domain D1.The resolution jump from the grid spacing of ERA-Interim data (0.75 • × 0.75 • ) to the one of the external WRF grid (12 km × 12 km) ranges from about 5.5 : 1 to 7 : 1 (in longitude and latitude, respectively).Its use is justified by Beck et al. (2004) who concluded that the considered one-way nesting strategies are comparably successful for precipitation simulations even if a large resolution jump of 10 : 1 is employed.Similarly, Liu et al. (2012) showed that downscaling ratios of 7 : 1, 5 : 1, and 3 : 1 in an WRF model perform better than smaller ratios in terms of rainfall.Resolution jumps (from the forcing data to the external model grid) higher than 5 : 1 and up to 10 : 1 have been successfully employed in various numerical studies (e.g., Galanis et al., 2006;Louka et al., 2008;Katsafados et al., 2011;Stathopoulos et al., 2013).The selection of the above-mentioned data sets was made for homogenization purposes of initial and boundary conditions for all cases.The upper air fields have been made available at 1000, 925,850,700,500,400,300,250,200,150,100,70, and 50 hPa.The sea-surface temperatures (SSTs) were derived from the daily National Centers of Environmental Prediction (NCEP) SST files at a very high horizontal resolution of 0.083 The Ferrier (Ferrier et al., 2002), RRTMG (rapid radiative transfer model application for global climate models) (Iacono et al., 2008), the Monin-Obukhov (Eta) (Janjic, 1996) and the Mellor-Yamada-Janjic level 2.5 (Mellor and Yamada, 1982;Janjic, 2002) schemes were used in all nests to represent microphysics, long-wave/shortwave radiation, surface layer, and boundary layer, respectively.The soil processes are represented by the NOAH (National Centers for Environmental Prediction/Oregon State University/Air Force/Hydrologic Research Lab) Unified model (Chen and Dudhia, 2001) in four layers (0-10, 10-40, 40-100, and 100-200 cm).Cumulus convection was parameterized only in nests D1 and D2 by the most recent version of the Kain-Fritsch scheme (Kain, 2004).
The model was initialized at 00:00 UTC, on 17 November 2007 (∼ 21 h before the event) for the 2007 tornado event.
For the tornado events of 2010 and 2011, the model was initialized at 00:00 UTC on 12 February 2010 (∼ 17 h before the event) and at 00:00 UTC on 20 September 2011 (∼ 15 h before the event), respectively.The time of initialization was selected in order to examine the model predictability more than 12 forecast hours ahead of the event and represent the values of previous diagnostic variables without spin up problems.
The model outputs were available at a 10 min interval.This interval was chosen so as to investigate the ability of the WRF to simulate, in short time intervals, the evolution of convective weather events influenced by topography modifications.

Tornadic events and synoptic analyses
In this section, a brief synoptic analysis for each tornado event is presented.Synoptic analyses concern the isobaric level of 500 hPa, based on ECMWF ERA-Interim analysis archive data set, and the mean sea level (MSL) pressure based on UK Met Office (UKMO) MSL pressure analysis.Figure 2 illustrates a collection of images revealing the force, the impact, and damage caused by each tornado to the local population.Hereafter, for reasons of brevity, we coded the tornadoes as TXX, where T stands for tornado and XX represent the last two digits of the year when tornado occurred (e.g., T10 = Tornado event that occurred in 2010).
The first tornado event (hereafter T07) formed at approximately 21:20 UTC (±5 min) on 17 November 2007, over Loutoufi (lat.: 38.28 • N, long.: 23.28 • E; Fig. 1), a small village located 9 km SW of the city of Thebes.The tornado dissipated in the NW Thebes urban area (Piri district), its path was 10 km, and scattered broken branches of olive trees were dispersed all over the area (Fig. 2a and b).Additionally, significant structural damage was documented in the village of Loutoufi and city of Thebes.Regarding the tornado force scale, it could be characterized as a T4-T5 on the TORRO scale, based on a tornado damage survey.Synoptic analysis at the isobaric level of 500 hPa at 18:00 UTC (not shown), depicts a closed cyclonic circulation over central Italy, implying a SW upper air stream over west Greece.A long trough line is positioned from southern Italy to northern Libya, accompanied by a thermal trough of −20 • C. UKMO MSL pressure analysis at 18:00 UTC (Fig. 3a) illustrates a cold front activity and an instability line over west Greece.
The second tornado event (hereafter as T10) occurred 2.5 km south of Chalkidiki, at the village of Vrastama (lat.: 40.36, long.: 23.54; Fig. 1), a non urban area 45 km south- east of Thessaloniki in northern Greece, on 12 February 2010.The T10 developed approximately between 17:10 and 17:35 UTC, and caused significant damage to a nearby greenolive processing unit and to several olive-oil-producing farms (more than 100 olive trees were uprooted).Based on the damage caused (Fig. 2c and d), it could be characterized as T4-T5 on the TORRO scale.An ECMWF ERA-Interim upper air analysis at the isobaric level of 500 hPa (not shown), during the day of T10 event, showed a closed cyclonic circulation associated with cold air masses (−35 • C), inducing a SW upper air flow over the area of interest at 12:00 UTC.The cyclonic circulation extended throughout the lower troposphere; UKMO MSL pressure analysis at 18:00 UTC (Fig. 3b) showed a cold front extended from northern Greece to Peloponnisos (southern Greece) and northern Libya.
The last tornado event (hereafter as T11) developed at approximately 15:20 UTC (±10 min), on 20 September 2011, over Vlychos (lat.: 38.68 • N, long.: 20.69 • E; Fig. 1), a small village on the island of Lefkada, over western Greece.T11 caused significant damage in Vlychos' marina (Fig. 2e and f) and one casualty was recorded.Tornado intensity has similarly been categorized as T4-T5 on the TORRO scale.At 500 hPa, a closed cyclonic circulation over central Italy (not shown) accompanied by a long wave trough along southern Italy, caused a SW upper air stream over western Greece.Mean sea level pressure analysis at 12:00 UTC illustrates frontal wave activity along western Greece (Fig. 3c).A detailed synoptic analysis accompanied by remote sensing data (satellite and radar images) for T10 event was also presented by Matsangouras et al. (2011a).The synoptic weather conditions and analysis during the day of the T07 event were similarly discussed by Matsangouras et al. (2012).Visualizations of weather conditions (spatial distribution of lightning activity and ground observations) for T10 and T11 events are presented in Fig. 4, based on the Hellenic National Meteorological Service (HNMS) archive data set.Images were produced via HNMS's METshell software (version v.2010.a)tailored for operational forecasting usage (Basiakos et al., 2000).

Model verification
The analysis of model errors in the control simulations performed on the day of each tornado is depicted in Table 2.The spin-up time (first 6 h) is not taken into account.The statistics are calculated using the available meteorological aerodrome reports (METARs) (from the meteorological stations of HNMS that are nearest to each tornado event) and the corresponding model values of the innermost domain (D03).The model output fields of the grid point closest to each station are used in the verification.Figure 1 illustrates the location of tornadoes and the meteorological stations (red and blue text, respectively).The distance of these stations from the corresponding tornado events range from 26 to 65 km and are considered representative of the meteorological conditions observed in the environment of tornadoes.Statistical measures of 10 m wind speed (WS), 2 m temperature (T ), 2 m dew point (Td), and mean sea-level pressure (MSLP) at the meteorological stations are depicted in Table 2.
The model overestimated the 10 m wind speed across all stations from 0.68 to 2.36 m s −1 .The MAE and RMSE of the 10 m wind speed ranged from 1.83 to 3.37 m s −1 and from 2.50 to 4.01 m s −1 , respectively.The 2 m temperature and the 2 m dew point were overestimated in all control simulations.The MAE (RMSE) of temperature was between 1.71 K (2.22 K) and 5.47 K (6.11 K) with the highest errors appearing at Tanagra airport (LGTG) (case of 17 November 2007).The dew point MAE and RMSE ranged from 1.19 to 3.37 K and from 1.60 to 4.42 K, respectively.Finally, the bias of the mean sea-level pressure was either close to 0 or negative, and its values fluctuated between 0.28 and −2.49hPa.Its MAE was between 1.26 and 2.49 hPa, while the range of the RMSE was from 1.43 to 2.64 hPa.
The above-mentioned values are in general agreement with statistical evaluations of modern high-resolution numerical weather prediction models in Greece.Gofa et al. (2008) validated the operational HNMS forecasts of SKIRON (0.06 • × 0.06 • ) and COSMO-GR (0.0625 • × 0.0625 • ) models for 2007 and found RMSE between: (a) 2.2 and 2.5 m s −1 for the 10 m wind speed, (b) 1.9 and 2.6 K for 2 m temperature, (c) 2.3 and 3.5 K for 2 m dew point, and (d) 1 and In this study, the model errors are generally within the values that appear in the literature, with a few exceptions (e.g., the 2 m temperature at LGTG and the 2 m dew point at LGPZ).The only parameter that systematically exhibited high values of RMSE is the 10 m wind speed.Some possible reasons for these errors are as follows: a.They correspond to high-impact weather events at individual stations and not to different weather types for multiple stations and long periods.In this study, most of the stations were under the influence of cumulonimbus activity for several hours.Even if a numerical model is successful in representing convective activity, it is very difficult to adequately represent the mesoscale flow induced by the storm.Ebert (2008) argues that although high-resolution numerical weather predictions can look quite realistic and be very useful, they have a "difficulty of predicting an exact match to the observations".
b.The model values were extracted from domain 3 that is integrated with a time step of 8 s.On the contrary, the observed 10 m wind speeds are 10 min average values.
c.The evaluation is based on METAR and not on SYNOP (surface synoptic observations).Thus, maximum errors up to +0.9 hPa in mean sea-level pressure and ±0.5 K in 2 m temperature/2 m dew point can be introduced for the statistical values because of the truncation and rounding  of the actual measurements of the above parameters, respectively, performed with METAR reports.Moreover, the fact that the METAR reports are available in hourly and half-hourly intervals, contrary to the 6-hourly intervals of SYNOP, which were employed in the abovementioned studies, allows the identification of temporary model errors, which could have been missed if the model validation was based on the latter reports.

Absolute vorticity budget
The vorticity field and the vorticity equation are useful tools, commonly employed in studies of atmospheric dynamics.
The absolute vorticity and the vorticity equation terms were derived at 10 min intervals at 1 km a.s.l.This height was selected because in nature it is generally located near the base of thunderstorms and in the simulations it was detected above the model topography (at the location of each actual tornado).
The values of the above terms, but not the conclusions, were modified when the height moderately increased (up to 2 km).The calculations were performed using the Grid Analysis and Display System (GrADS) and were based on the simulated WRF parameters of the innermost domain (D3).The necessary upper air fields were vertically interpolated at 250 m intervals before the calculation of the vertical derivatives.All the derivatives were estimated using second-order centered finite differences.
Figures 5-7 illustrate time series of the maximum values of the absolute vorticity and the vorticity equation terms of horizontal and vertical advection, tilting/twisting, and divergence, in a box of 0.2 • × 0.2 • latitude-longitude centered at the actual location of each tornado.It was chosen to exhibit the results in a limited area and not locally, in order to consider any shift of the meteorological conditions due to model error.Both experiments (with and without topography) are shown for each tornado case.The solenoidal term is not displayed because it was computed to be at least four orders of magnitude less than the other terms.A very limited number of grid points with topography above or equal to 1 km within the above-mentioned area (0.2 • × 0.2 • ) were not taken into account.
The most important term of the vorticity equation appears to be the horizontal advection, with maximum values up to 0.35 (0.29) s −1 h −1 in the case of 12 February 2010 with (without) topography.The predominance of this term in these case studies is likely to be associated with the prevailing strong synoptic forcing due to the cold front.The vertical advection and twisting/tilting of horizontal vorticity to the vertical direction appear to be similarly important without significant differences, while the divergence term has generally a smaller contribution than the previous two terms.In real supercells, the tilting of horizontal vorticity generated along the gust front into vertical vorticity is a significant source of low-level vorticity (Houze, 1993).However, this mechanism cannot be adequately represented in the simulations of this article, despite their fine (but not fine enough) grid spacing (1.333 km × 1.333 km in D3).In all cases, the simulations with topography are associated with higher values of maximum absolute vorticity in the area of interest than those without topography (Figs.5-7).Similarly, the vorticity equation terms that determine the Eulerian time change of vorticity display higher maximum values in the experiment with topography.The differences between the two experiments, and the concomitant importance of topography for tornadogenesis, are more distinguishable in the cases of 17 November 2007 (Fig. 5) and 12 February 2010 (Fig. 6).The former tornado occurred at a location surrounded by mountains and the latter one developed in the mountains of Chalkidiki in northern Greece.

Results and discussion
In this section, we present the results of numerical simulation analyses of EHI, BRN, SRH, and MCAPE diagnostic variables based on the actual topography and modified topogra-  phy in the inner domain box (D3 as described in Sect.2).The actual topography is illustrated in meters in Fig. 1a.Hereafter, for brevity reasons, the numerical simulation results based on modified (reducing) topography of D3 domain will be coded as TOPOMX, where X stands for the percentage (%) of modification.Thus, TOPOM0 means modified topography by 0 % (actual topography), and TOPOM100 refers to modified topography by −100 % (without topography).Figure 8a-d, illustrate TOPOM0 (solid lines) and TOPOM100 (dotted lines) simulation results 90 min prior to and 60 min after tornadogenesis for EHI, SRH, BRN, and MCAPE across every tornadic case study.All values for the two sets of simulations (TOPOM0 and TOPOM100) depict the values of diagnostic indices in radii less than 20 km from the tornadogenesis location.The spatial distribution of the four diagnostic variables of TOPOM0 simulations at the time of tornadogenesis is presented in Fig. 9 for each tornado case.
The EHI analysis of T07 TOPOM0 numerical simulation reveals significant high values (EHI > 6) prior to tornadogenesis (60 min prior) that gradually increase (EHI = 8) at the time of tornadogenesis (Fig. 8a, solid green line with squares).Figure 9a  less than 15 km to the south), while 30 min after tornadogenesis EHI reduces to 6 (Fig. 8a).In contrast to TOPOM0, TOPOM100 (Fig. 8a, dotted green line with squares) EHI analysis depicts a gradually reduction of EHI values (from 5 to 3) prior to the time of tornadogenesis, documenting the influence of topography.During T10 numerical simulations with TOPOM0 (Fig. 8a, solid red line with triangles) and TOPOM100 (Fig. 8a, dotted red line with triangles), EHI analysis presents significant high values (EHI > 5) for both simulation sets.More specifically, the TOPOM0 experiment shows a gradual increase of the EHI value (from 5 to 8, 80 min prior to tornadogenesis) (Fig. 8a). Figure 9b depicts this value too close to the tornadogenesis location (less than 10 km to the WSW), while after tornado formation the EHI declines to 5. The EHI analysis of T10 TOPOM100 experiment presents lower values (compared to TOPOM0) that started reducing prior to tornadogenesis.Regarding the T11 (a coastal tornado event) TOPOM0 (Fig. 8a, solid blue line with circles) and TOPOM100 (Fig. 8a, dotted blue line with circles) numerical simulations, EHI analysis remains identical, suggesting that topography did not influenced EHI.
The spatial distribution of EHI in T11 TOPOM0 experiment shows a value of EHI > 2, 30 km SE of tornado location (Fig. 9c).The SRH analysis of T07 TOPOM0 (Fig. 8b, solid green line with squares) numerical simulation illustrates that maximum values fluctuate between 1000 and 1200 m 2 s −2 prior to and after tornadogenesis.Due to the eastern propagation of a cell area with high values, this is depicted in Fig. 9e as less than 20 km SE from the tornado location.TOPOM100 (Fig. 8b, dotted green line with squares) numerical simulation spatial distribution depicts the same cell area with a value of 200 m 2 s −2 less than SRH TOPOM0, suggesting an influence of topography at this tornado study.The influence of topography is also evident in the SRH analysis of T10 case study.T10 TOPOM0 (Fig. 8b, solid red line with triangles) numerical simulation suggests that the maximum values were constant at 1200 m 2 s −2 (prior to tornadogenesis) and increased significantly (2000 m 2 s −2 ) as the storm propagated eastwards.Figure 9e illustrates the maximum SRH T10 TOPOM0 value (1200 m 2 s −2 ) in less than 10 km WSW.T10 TOPOM100 (Fig. 8b, dotted red line with triangles) results follow the TOPOM0 spatial distribution pattern (significant increase of SRH after tornadogenesis) but in lower values (200 m 2 s −2 less).Regarding T11 simulations, SRH analysis does not reveal any influence of topography, but dis-plays a similar output to the EHI analysis.From these two outputs, regarding T11 simulations, we become aware that similar CAPE numerical investigation is not going to depict any positive or negative influence of topography as EHI is expressed in terms of SRH and CAPE (Table 2).
Regarding BRN analysis of T07 TOPOM0 (Fig. 8c, solid green line with squares), the maximum values of diagnostic variables show an incremental trend (from 400 to 600 m 2 s −2 ) prior to tornadogenesis and a declining one after tornadogenesis (from 600 to 200 m 2 s −2 ).The spatial distribution of BRN TOPOM0 simulation at the time of tornadogenesis illustrates an area of maximum values (600 m 2 s −2 ) west (< 10 km) of the tornado location (Fig. 9g).Concerning TOPOM100 (Fig. 8c, dotted green line with squares) the BRN analysis follows the TOPOM0 pattern (increase of values prior to tornadogenesis and decrease of values later) but the maximum values fluctuate between 200 and 400 m 2 s −2 .
The latter BRN analysis for T07 reveals the influence of topography as the TOPOM100 simulation demonstrates lower values compared with the TOPOM0 simulation.Studying the two sets of T10 numerical simulations, TOPOM0 (Fig. 8c, solid red line with triangles) and TOPOM100 (Fig. 8c, dotted red line with triangles), we come to the same conclusion regarding the analysis of BRN diagnostic variable.The TOPOM0 BRN spatial distribution illustrates higher values (Fig. 9h, 300 m 2 s −2 , 10 km west of the tornado location) than TOPOM100.In contrast to T10, the two numerical experiments of T11 do not reveal any influence of topography over BRN analysis.In both simulations of T11 (Fig. 8c) BRN variable analysis displays an identical spatial distribution.
At T07 TOPOM0 simulation, MCAPE analysis (Fig. 8d, solid green line with triangles) shows that values remain steadily higher than 1500 J kg −1 .Regarding T07 TOPOM100 simulation, MCAPE analysis (Fig. 8d, solid green line with triangles) illustrates a gradual decrease of MCAPE prior to tornadogenesis.MCAPE spatial distribution (Fig. 9k) from TOPOM0 simulation at the time of tornadogenesis illustrates an area of high values south of tornado location.Regarding T10, MCAPE analysis of TOPOM0 (Fig. 8d, dot red line with triangles) and TOPOM100 (Fig. 8d, dotted red line with triangles) simulations indicates that topography does not significantly affect the evolution of this diagnostic variable.The MCAPE at TOPOM100 analysis is reduced only for a short time frame (from 30 min prior to 20 min after tornadogenesis).MCAPE analysis for the simulation sets of T11 illustrates an identical pattern of maximum values, depicting that topography does not influence the MCAPE diagnostic variable in this tornado event.
Based on EHI, SRH, BRN, and MCAPE analyses of TOPOM100 and TOPOM0 simulations, the complex terrain was important in T07 and T10, with the former exhibiting higher values.Topography appeared to constitute an unimportant factor during T11 event, utilizing EHI, SRH, BRN, and MCAPE diagnostic variables.Conversely in the T07 event, topography was an important factor across all diagnostic variables, since EHI and MCAPE variables at TOPOM0 simulation actually doubled compared to the TOPOM100 simulation at the time of tornadogenesis.During the T07 event, the storm propagated ENE, over a steep slope terrain, from the coast at the Gulf of Corinth to a complex terrain of 1200 m elevation.EHI, MCAPE and SRH variables showed a significant increase 30 min prior to the T07 tornado event.

Summary and conclusions
Three tornado events that occurred in recent years in Greece were selected for numerical experiments, in order to investigate the role of topography in significant tornadogenesis events triggered under strong synoptic scale forcing.The first tornado event (T07) affected Thebes (Boeotia, on 17 November 2007), the second event (T10) Vrastema (Chalkidiki, on 12 February 2010) and the last one (T11) Vlychos (Lefkada, on 20 September 2011).These events have been associated with synoptic scale forcing, while their intensities were T4-T5 (TORRO scale) causing significant damage.
Considering the time of tornado occurrence and MSL synoptic analysis charts, from UK Met Office (UKMO), all tornadoes were formed within storms that were associated with a synoptic force ingredient.All events occurred in the profrontal activity area, ahead of cold front.
Numerical simulations were performed using the nonhydrostatic weather research and forecasting (WRF, ARW core) model.Three one-way nested domains were utilized with spatial resolution of 12 km for D1, 4 km for D2 and 1.333 km for D3 domain.Initial and lateral boundary conditions for the D1 domain were obtained from ECMWF ERA-Interim data set, with telescoping nested grids that allow for the representation of atmospheric circulations ranging from the synoptic scale down to the mesoscale.The WRF-ARW model setup was able to forecast significant values of diagnostic variables with a lead time of more than 12 h.In all experiments the topography of the inner grid (D3) was modified by: (a) 0 % (actual topography) and (b) −100 % (without topography).Analyses concerned a data set of diagnostic instability variables: EHI, BRN, SRH, and MCAPE.
The model verification performed in this study has shown that model errors are generally within the values that appear in the literature, with a few exceptions (e.g., the 2 m temperature at LGTG and the 2 m dew point at LGPZ).The analysis of absolute vorticity budget revealed that in all cases the simulations with topography are associated with higher values of maximum absolute vorticity in the area of interest than those without topography.Similarly, the vorticity equation terms that determine the Eulerian time change of vorticity display higher maximum values in the experiment with topography.The differences between the two experiments, and the concomitant importance of topography for tornadogenesis, are more distinguishable in the cases of 17 November 2007 and 12 February 2010.The similar values of maximum absolute vorticity and maximum vorticity equation terms in the two experiments of T11 event suggested that, in this case, topography played a minor role, while the synoptic forcing dominated.
Numerical simulations revealed that complex topography constituted an important factor during 17 November 2007, 12 February 2010 events, based on EHI, SRH, BRN, and MCAPE analyses.Conversely, topography around the 20 September 2011 event was a less important factor based on EHI, SRH, BRN, and MCAPE analyses.
Further investigation of additional cases or idealized simulations is suggested in order to delineate the role of topography in tornadogenesis.Our work towards this aim is currently underway.

Figure 1 .
Figure 1.The upper image (a) illustrates the actual topography of Greece with the location (black dots) of tornadoes and meteorological stations (red and blue abbreviations, respectively).The three nested domains used for numerical investigations of tornado events (lower images) are illustrated during (b) 17 November 2007, (c) 12 February 2010 and (d) 20 September 2011.In (b), (c) and (d) images D1 is the outer domain, D3 is the inner domain and D2 is the intermediate domain.

Figure 3 .
Figure 3. Mean sea level analysis charts, from the UK Met Office, depicting the frontal activity around the tornadoes' time of occurrence for T07 (a), T10 (b) and T11 (c) tornado events.In every image, red triangles indicate the tornado locations.

Figure 4 .
Figure 4. Visualization of lightning distribution in a 15 min step, associated with the T11 (upper images) and T10 (lower images) events.Red and yellow stars represent the spatial distribution of lightning activity 15 and 30 min prior tornado formation, respectively.Tornadoes developed within the above hourly time window and the red X in a red circle denotes the location.METAR data from the closest weather stations, Preveza (LGPZ) and Thessaloniki (LGTS) for T11 and T10, respectively, were also presented.

Figure 5 .
Figure 5.Time series of WRF maximum simulated values of Absolute Vorticity (s −1 ) and of the vorticity equation terms (s −1 h −1 ) of horizontal (Hadv) and vertical advection (Vadv), tilting/twisting (TT), and divergence (DIV), at 1 km a.s.l. in a box of 0.2 • × 0.2 • latitude-longitude centered at the actual location of the tornado of 17 November 2007, (a) with topography and (b) without topography.Time is in UTC.

)Figure 6 .
Figure 6.Similar to Fig. 5 but for the case of 12 February 2010.

Figure 7 .
Figure 7. Similar to Fig. 5 but for the case of 20 September 2011.

Figure 8 .
Figure 8. Graphics (a-d) illustrating TOPOM0 (solid lines) and TOPOM100 (dotted lines) simulations results 90 min prior to and 60 min after tornadogenesis for EHI (a), SRH (b), BRN (c), and MCAPE (d) for every tornadic case study.The results of T07, T10, and T11 simulations are illustrated with green lines with squares, red lines with triangles, and blue lines with circles, respectively.All values for two sets of simulations (TOPOM0 and TOPOM100) illustrate the values of indices in radii of less than 20 km from the tornadogenesis location.

Figure 9 .
Figure 9. Spatial distribution of the four diagnostic variables of TOPOM0 simulations at the time of tornadogenesis for each tornado case (T07 left column, T10 middle column and T11 right column.

Table 1 .
A selection of indices commonly used for severe storm forecasting.In the formulae, T denotes temperature and D denotes dew point in • C, with a subscript indicating at what mandatory pressure level

Table 2 .
Statistical measures of 10 m wind speed (WS), 2 m temperature (T ), 2 m dew point (Td) and mean sea-level pressure (MSLP) at the meteorological stations of the Hellenic National Meteorological Service at the airports of Tanagra (LGTG), Thessaloniki (LGTS), Cephalonia (LGKF), and Preveza (LGPZ).The model output fields of the D03 grid point closest to each station are used.ME = mean error (WRF-simulated minus observed value), MAE = mean absolute error, RMSE = root mean squared error, sample size = number of valid pairs of simulated and observed values.The date of the model validation and the distance of each station from the location of the corresponding tornado event are also indicated.The spin-up time (first 6 h) is not taken into account.