Satellite-Based Study and Numerical Forecasting of Two Tornado Outbreaks in the Ural Region in June 2017

Strong tornadoes are common for the European part of Russia but happen rather rare east of the Urals. June 2017 became an exceptional month when two tornado outbreaks occurred in the Ural region of Russia, yielded $3 million damage, and resulted in 1 fatality and 14 injuries. In this study, we performed detailed analysis of these outbreaks with different data. Tornadoes and tornado-related environments were diagnosed with news and eyewitness reports, ground-based meteorological observations, sounding data, global numerical weather prediction (NWP) models data, synoptic charts, satellite images, and data of specially conducted aerial imaging. We also estimated the accuracy of short-term forecasting of outbreaks with the WRF-ARW mesoscale atmospheric model, which was run in convection-permitting mode. We determined the formation of 28 tornadoes during the first outbreak (3 June 2017) and 9 tornadoes during the second outbreak (18 June 2017). We estimated their intensity using three different approaches and confirmed that, based on the International Fujita scale (IF), one of the tornadoes had the IF4 intensity, being the first IF4 tornado in Russia in the 21st century and the first-ever IF4 tornado reported beyond the Ural Mountains. The synoptic-scale analysis revealed the similarity of two outbreaks, which both formed near the polar front in the warm part of deepening southern cyclones. Such synoptic conditions yield mostly weak tornadoes in European Russia; however, our analysis indicates that these conditions are likely favorable for strong tornadoes over the Ural region. Meso-scale analysis indicates that the environments were favorable for tornado formation in both cases, and most severe-weather indicators exceeded their critical values. Our analysis demonstrates that for the Ural region, like for other regions of the world, combined use of the global NWP model outputs indicating high values of severe-weather indices and the WRF model forecast outputs explicitly simulating tornadic storm formation could be used to predict the high probability of strong tornado formation. For both analyzed events, the availability of such tornado warning forecast could help local authorities to take early actions on population protection.


Introduction
Tornadoes are one of the most devastating and hard-to-predict atmospheric phenomena observed in mid-latitudes. Moreover, the threat of tornadoes in Europe and Russia is substantially underestimated [1][2][3]. In particular, it was considered for a long time that tornadoes occur mainly in European Russia and are extremely rare east of the Urals [4]. Recently, Chernokulsky et al. [2] showed that tornadoes are common for both European and Asian parts of Russia but confirmed that strong tornadoes (≥F3 on the Fujita scale [5]) are substantially less frequent east of the Urals. Particularly, We examined synoptic-scale and mesoscale environments of two tornado outbreaks' formation using different data, including numerical experiments with the WRF model. We estimated the model performance to forecast the events as well as model sensitivity to grid size and forecast lead time (from 12 to 36 h). Since Doppler radar data were unavailable for the study area, the validation of WRF simulation results was performed with the Meteosat-8 satellite images and satellite-based estimates of tornado-induced disturbances in forests. We describe the experimental design in Section 2. Section 3 presents the results of synoptic-scale and mesoscale analysis of the events and analysis of the tornadoes' impact, while Section 4 provides the main results of the numerical simulation of tornadic storms. Section 5 summarizes the main findings of the paper.

Overview of the Data Sources
We used various data sources to obtain the most complete information on the characteristics of the analyzed severe weather outbreaks of 3 and 18 June 2017. These data are summarized below.
Ground-based observations included regular observations from the weather stations of Russian state weather service (69 stations in the Sverdlovsk, Tyumen, and Kurgan regions) and sounding data obtained from the University of Wyoming database [50] that were selected based on the proximity sounding approach [51]. For the event of 3 June 2017, we used the 1200 UTC sounding of Verhnee Dubrovo station (WMO ID: 28445). However, we failed to find suitable sounding data for the event of 18 June 2017.
Severe weather reports included compiled eye-witness photo, video, and damage reports from the media, social networks, and from the European Severe Weather Database [49] and were used to clarify the tornadoes' impact, path, and intensity.
Satellite data included medium-resolution images from the Landsat-8 and Sentinel-2 satellites and high-resolution images from Google Earth, which were both used to identify forest damage tracks induced by windstorms and tornadoes, and estimate the path length, width, and intensity of the tornadoes (see Section 2.2 for more details). Additionally, we used the Meteosat-8/SEVIRI images, obtained from the EUMETSAT Earth Observation portal [52], to explore convective storms' evolution, identify the signatures of intense updrafts, and validate the WRF model skills in storm track forecasting. We used the RGB (Red-Green-Blue) combination of the High-Resolution Visible (HRV) and 10.8 µm infrared (IR10.8) spectral bands, known as "HRV cloud", and also cloud top temperature estimated by the IR10.8 band, to detect overshooting tops (OTs), cold-ring, and cold U/V (CRCUV)-shaped features on the cloud top. According to previous studies [12,[53][54][55][56], OTs and CRCUV-shaped signatures indicate strong updrafts and mesocyclones in mesoscale convective systems (MCSs).
In addition, we obtained high-resolution aerial orthophotographs (with 0.2 m pixel size) of the damage track of the tornado that passed near Visim settlement on 3 June 2017. Aerial imagery was performed using a quadcopter on 18 May 2018 along the tornado track. In total, an area of 20 km long and 1 km wide was covered by orthophotographs.
Synoptic-scale fields of atmospheric variables were obtained from the short-term forecast data of the GFS (Global Forecast System [57]) and GEM (Global Environmental Multiscale Model [58]) NWP models. Synoptic charts obtained from the Roshydromet archive were also utilized. The GFS and GEM models forecasts (with~25 km grid size) from 0000 UTC 03 June 2017 and 0000 UTC 18 June 2017 were downloaded from the NCEP [57] and CMC [58] ftp-servers, respectively.
The WRF mesoscale NWP model with the ARW (Advanced Research WRF) dynamic core version 3.8.1 (WRF-ARW) [39] was used for explicit simulation of MCSs (see Section 2.3 for details).

Tornado Tracks Identification and Characteristics Estimate
Initial data on the studied tornado events was obtained from eye-witness observations and damage reports, presented in news and social media and then further uploaded to the European Severe Weather Database [49] by one of the authors (IA) of this study. However, the intensity of tornadoes was estimated preliminarily, while the characteristics of the tornado paths were mostly undetermined. We examined them in this paper using various satellite observations. In particular, we used satellite data in three-stage searching of tornado damage tracks in forests of storm-affected regions.
For the first stage, we used pre-and post-event Landsat-8 and Sentinel-2 images downloaded from the USGS Data Archive [61]. We applied the previously proposed method of tornado track identification in forested regions [10,62,63] and delineated numerous windthrow areas caused by tornadoes and storm winds. Additionally, we used the Global Forest Change data [64] to confirm the identified forest loss areas.
For the second stage, we analyzed high-resolution satellite images available on Google Earth to confirm or reject potential tornadic nature of windthrow areas, based on its typical signatures [13,62,63]. Most of the storm-affected region was covered by summertime high-resolution images in 2017-2018 (after the storm events), which allowed to determine the cause of windthrow, i.e., discriminate between a tornado and squall. Using high-resolution images, we visually inspected in detail areas adjacent to storm-and tornado-damage tracks identified by Sentinel-2 images and found several additional tornado tracks. Some of them partially overlapped with squall-induced windthrow areas.
For the third stage, we used the Meteosat-8/SEVIRI images to improve the searching of tornado tracks, which could be missed after the first and second stage. We searched OTs and CRCUV-shaped features on the cloud tops, which indicate intense updrafts [12,[53][54][55][56], and delineated their tracks. Along these tracks, we performed additional visual inspection of local windthrow areas using high-resolution images from Google Earth. This procedure helped us to discover several previously unknown tornado tracks.
We used high-resolution images to determine tornado/squall path geometrical characteristics. In the absence of high-resolution images, we used the Sentinel-2 data. We determined the path length, maximum width, and damaged area. More details on the procedure of geometrical characteristics determination can be found in [62,63].
For tornado events, we estimated tornado intensity along with geometrical characteristics. Most of the tornadoes passed outside of settlements and were found only with the use of satellite images of forest damage tracks. Ideally, forest damage survey or high-resolution aerial imagery of damaged forest is needed to correctly estimate the intensity of such tornadoes [13,15,65,66]. Due to their lack, we used high-resolution (with~0.5 m pixel size) satellite images from Google Earth to assess tornado intensity. The spatial resolution of satellite images is several times lower than the same of aerial images; however, they can provide valuable information to estimate tornado damage. In particular, we used three different approaches.
The first approach of tornado intensity estimation is based on the recommendations for the International Fujita (IF) Scale [65], which were elaborated by the European Severe Storm laboratory, and can be applied to Russian forest species and, with some assumptions, to Russian buildings. It was primarily applied for tornadoes covered by photo or video damage reports. With substantial limitations, this method may be applied to tornadoes identified using high-resolution satellite images of damaged forests. Particularly, we examined the forest damage degree that may be of partial or total canopy removal (less and more than 75% of windthrow area), which corresponds to the IF0 and ≥IF1 tornado Atmosphere 2020, 11, 1146 5 of 34 intensity, respectively [65]. It is important to note that all widespread tree species in tornado-damaged forests (spruce, pine, birch, and aspen) are classified as low firmness species (softwood) according to [65]. Therefore, tornado having intensity ≥IF1 is able to destroy forest stands. However, we did not able to determine whether the trees were uprooted, snapped, or debarked because of the tornado, which can indicate higher intensity levels [65]. The one exception is the tornado that passed near Visim and was covered by high-resolution aerial images. We found this the tornado caused total canopy removal and snapped most coniferous trees within the damage track, while deciduous trees were mostly felled down; as a result, we rated it as IF2 [65].
The second method of the enhanced Fujita (EF) scale intensity determination is based on the approach proposed by Peterson et al. [66]. Within this approach, the EF-scale level is a function of forest damage severity, which is expressed as the percentage of uprooted or snapped trees on 100 × 100 m areas. For instance, the EF0 level corresponds to 18-38% of areas of blown trees (10-48% within the 95% confidence interval) (see Figure 7 from [66] for details). We applied this approach for tornadoes that damaged forest stands and for which we had post-event high-resolution satellite images. Note, we did not estimate the EF-rating based on recommendations [67] since they were proposed for U.S. construction practices. The approach [66] should also be used with caution since it was developed using linkages between forest wind resistance and wind speed that were obtained for forest species in North America rather than northern Eurasia.
The third approach allows to determine minimum probable Fujita-scale intensity based on tornado path length and width [2,10,62]. We used Weibull distribution parameters from [68] that link tornado intensity with path length and maximum width. The approach is described in full in [62] (see particularly their Figure 6). It was applied for all tornadoes whose path characteristics were determined.
Combining three approaches, we assigned the most probable intensity rating (on the IF-scale) to all tornadoes that occurred on 3 June and 18 June 2017 that were restored using the aforementioned approaches. The Appendix A presents the list of these tornadoes (Table A2) and also the list of strong wind events reported at weather stations and based on forest damage analysis (Table A3).

WRF Model Configuration and Forecast Accuracy Assessment
We used the WRF-ARW v.3.8.1 model to simulate the analyzed events. This is a non-hydrostatic model with several parameterization options [39,69]. We performed several numerical experiments with different horizontal grid spacing, with and without a nested grid. We used the GFS-NCEP model forecast with a 0.25 • grid size as initial and boundary conditions. The model was run in a convection-permitting mode, which was previously found to be optimal for tornadic storm simulation [29,33,41,42,45,46]. We used explicit convection-permitting simulation since most of the simulated convective storms were meso-α-scale events. The nudging procedure was not implemented. Table 1 summarizes the WRF model settings, and Figure 1 shows the model domains.
The sensitivity numerical experiments with the WRF-ARW model included two main steps. On the first step, we estimated the influence of horizontal resolution on the results of tornadic storm simulation. The simulation was performed with the horizontal grid spacing of 9, 7.2, and 3 km. For the 9-km grid spacing, one nested domain with a 3-km resolution grid performed with the one-way nesting technique was included. The model assimilated the GFS-NCEP data for 0000 UTC of 3rd and 18th June 2017 for the first and second events and operated for 24 h till 0000 UTC of 4th and 18th June 2017, respectively. Since both tornadoes were observed around 1200 UTC, we could consider this experiment as a forecast with an approximately 12-h lead time. As a result, we determined the most appropriate model configuration (in terms of resolution) to simulate tornadic storms with this lead time (see Section 4 for more details).

Storm Events Description
In May and June 2017, a deep and unusually long-lived tropospheric trough dominated over the European Russia and caused a long-term period of cold and wet weather with average temperature below the norm by 2-4 °C. Arctic cold air masses spread to the Black, Caspian, and Aral Seas, inducing the formation of several southern cyclones that moved to the Ural region. These cyclones yielded favorable conditions for convective storm development, specifically on 3 June 2017 and 18 June 2017. These days, severe weather outbreaks with heavy rainfalls, strong wind ≥25 m s −1 , large hail up to 40 mm, and tornadoes formed in the Sverdlovsk and Kurgan/Tyumen regions, respectively. Both events caused substantial damage to settlements, agricultural crops, and forests. One fatality and at least 14 injuries were reported.

Large-Scale Synoptic Features
On 3 June 2017, the Ural region was influenced by the forward side of a deep tropospheric cyclone with a center over north of European Russia. The frontal zone stretched over the Ural Mountains from south to north with a temperature gradient up to 12 °C/500 km at 850 hPa between rather cold European Russia (−5 °C) and warm Western Siberia (up to 18 °C). A mid-tropospheric jet stream with wind speeds up to 35 m s −1 at 500 hPa stretched from the southern part of European Russia to Western Ural.
In the afternoon of 2 June 2017, a surface low (with SLP ~1000 hPa) formed on the polar front between Caspian and Aral Seas. After 0000 UTC 3 June, it began to move with speed around 16-22 m s −1 along the jet stream. At 1200 UTC 3 June, the low deepened (SLP amounted to 993 hPa) and centered near 55° N, 60° E, where it resulted in a severe weather outbreak ( Figure 2a). Afterward, the low moved northward to the Northern Ural and continued to deepen (up to 987 hPa at 0000 UTC 4 June 2017).
In the afternoon of 3 June 2017, the polar front stretched along the Ural Mountains. It divided cold air mass to the west (i.e., in the Perm region), with the maximum air temperature at 2 m AGL (Tmax) around +16°…+20 °C, and warm air mass to the east (i.e., in the Sverdlovsk and Chelyabinsk regions) with Tmax up to +31 °C. Severe convective storms with squalls, large hail, and tornadoes formed in the warm sector and on the cold front near the occlusion point. The synoptic situation of On the second step, we estimated the sensitivity of forecast accuracy to the lead time. The WRF model assimilated the GFS-NCEP data for 0000 and 1200 UTC of preceding days, i.e., 2 June 2017 and 17 June 2017, respectively. As a result, we obtained the forecasts of tornadic storms with 24-and 36-h lead times. For this experiment, we used the model configuration that provided a more accurate forecast on the first step. The Doppler radar system installed in Russia does not cover the region of interest yet [70]. Therefore, we used the Meteosat-8 SEVIRI satellite images to validate results of the WRF outputs. We compared satellite-observed and WRF-simulated characteristics of convective tornadic storms-path, timing, and cloud top temperature. Additionally, the location and timing of simulated storms was compared with the observed tornado tracks. Tornadic storms were identified based on sea-level pressure (SLP), 0-3 km SRH, maximum radar reflectivity, and 10 m wind speed fields. In particular, we considered that meso-scale low pressure systems (having 4-10 hPa depth), high values of SRH (>1000 m 2 s −2 ), and strong wind gusts (≥30 m s −1 ) with well-pronounced wind convergence indicated potential tornadic storms. Similar approaches used to identify convective storms based on WRF model outputs have been successfully used previously (see, e.g., [17,42,45]). We should stress that the use of SRH alone may not have meaningful application since simulated convective cells themselves significantly affect wind field and SRH. In theory, the updraft helicity may better indicate rotating updrafts (i.e., supercell mesocyclone) (see, e.g., [22,32,33]); however, due to some technical issues, we did not calculate this parameter in our analysis.

Storm Events Description
In May and June 2017, a deep and unusually long-lived tropospheric trough dominated over the European Russia and caused a long-term period of cold and wet weather with average temperature below the norm by 2-4 • C. Arctic cold air masses spread to the Black, Caspian, and Aral Seas, inducing the formation of several southern cyclones that moved to the Ural region. These cyclones yielded favorable conditions for convective storm development, specifically on 3 June 2017 and 18 June 2017. These days, severe weather outbreaks with heavy rainfalls, strong wind ≥25 m s −1 , large hail up to 40 mm, and tornadoes formed in the Sverdlovsk and Kurgan/Tyumen regions, respectively. Both events caused substantial damage to settlements, agricultural crops, and forests. One fatality and at least 14 injuries were reported.

Large-Scale Synoptic Features
On 3 June 2017, the Ural region was influenced by the forward side of a deep tropospheric cyclone with a center over north of European Russia. The frontal zone stretched over the Ural Mountains from south to north with a temperature gradient up to 12 • C/500 km at 850 hPa between rather cold European Russia (−5 • C) and warm Western Siberia (up to 18 • C). A mid-tropospheric jet stream with wind speeds up to 35 m s −1 at 500 hPa stretched from the southern part of European Russia to Western Ural.
In the afternoon of 2 June 2017, a surface low (with SLP~1000 hPa) formed on the polar front between Caspian and Aral Seas. After 0000 UTC 3 June, it began to move with speed around 16-22 m s −1 along the jet stream. At 1200 UTC 3 June, the low deepened (SLP amounted to 993 hPa) and centered near 55 • N, 60 • E, where it resulted in a severe weather outbreak ( Figure 2a). Afterward, the low moved northward to the Northern Ural and continued to deepen (up to 987 hPa at 0000 UTC 4 June 2017).
In the afternoon of 3 June 2017, the polar front stretched along the Ural Mountains. It divided cold air mass to the west (i.e., in the Perm region), with the maximum air temperature at 2 m AGL (Tmax) around +16 • . . . +20 • C, and warm air mass to the east (i.e., in the Sverdlovsk and Chelyabinsk regions) with Tmax up to +31 • C. Severe convective storms with squalls, large hail, and tornadoes formed in the warm sector and on the cold front near the occlusion point. The synoptic situation of the 3 June 2017 event generally corresponds to the second typical scheme of tornado formation over Northern Eurasia according to [4].
On 18 June 2017, the eastern part of the Ural region was influenced by a forward side of a tropospheric low centered over the Volga region. A waving polar front stretched from south to north over northern Kazakhstan and the Kurgan and Tyumen regions. The frontal zone in the middle troposphere was similarly oriented; wind speed at 500 hPa did not exceed 25 m s −1 . the 3 June 2017 event generally corresponds to the second typical scheme of tornado formation over Northern Eurasia according to [4].   (Figure 2b). In the next hours, its movement direction changed to northwestern and speed decreased to 8-11 m s −1 . At 0600 UTC 19 June, the low centered over the western part of the Tyumen region, and SLP decreased to 982 hPa.
Deep convection started near the occlusion point of the cyclone at 0930 UTC 18 June when the temperature gradient over the polar front was up to 12-14 • C/500 km at 850 hPa. As for the first event, the synoptic situation of this event corresponds to the second typical scheme of tornado formation over Northern Eurasia according to [4].

Diagnostic Variables, Estimated by the Global NWP Models, and Sounding Data
Several thermodynamic and kinematic variables (indices) related to the formation of tornadoes were calculated according to the GFS/NCEP and GEM/CMC data. We used the so-called proximity soundings approach [60] and found the maximum values of the severe-weather indices estimated within a 100-km area around the observed tornadoes. The selected indices are most frequently used to estimate the risk of tornadogenesis (see, e.g., [17,22,[40][41][42][43][44]59]).
According to the NWP-based calculating indices, on 3 June 2017 between 0900 and 1500 UTC, the environments were favorable for the development of severe thunderstorms and tornadoes in the Sverdlovsk and Chelyabinsk regions (Table 2, Figure 3). In particular, most of indices exceeded critical values. That is, CAPE ranged from 1500 to 2000 J kg −1 and low-level shear (LLS) reached 15 m s −1 according to the GEM data. These high values of LLS indicate the presence of a low-level jet stream that could contribute to the formation of tornadoes. Simultaneously, DLS was rather low for tornado occurrence (4 and 8 m s −1 according to the GFS and GEM, respectively). However, it is very likely that both NWP models substantially underestimated it, since DLS was 20 m s −1 according to the sounding data. Models showed the highest values of 0-3 km SRH (up to 480 m 2 s −2 ) near the axis of the trough oriented along the Ural Mountains. The maximum calculated EHI values exceeded 2.0 (up to 4.2 according to GFC), which is considered favorable for tornado occurrence [42]; STP values were also above the critical level (1.0); however, the SCP parameter lacked predictive power for this event ( Table 2). amounted to −72 J kg −1 . Thus, the atmospheric environments were exceptionally favorable for the development of a severe weather outbreak, including strong tornadoes.   The GFC-based soundings are sensitive to grid point selection. In particular, the sounding in the actual tornado point (Figure 3a) showed unusually high lifted condensation level. In fact, such conditions with dry air mass spreading from south-east in the low-layer were observed 100 km east of the tornado-affected area. It is of note that such environments resulted in dry squalls (severe wind gusts with almost no precipitation) in the central and southeastern parts of the Sverdlovsk region. The sounding data from the Verhnee Dubrovo station obtained at 1200 UTC 3 June 2017 confirmed the presence of low-layer dry air (CAPE was only 11 J kg −1 ) and indicated relatively high wind shear, with mid-tropospheric (at 511 hPa) wind speeds reaching 27 m s −1 , and a well-pronounced low-level, jet stream with wind speed reaching 18.5 m s −1 at the 879 hPa level. LLS and DLS amounted to 14 and 20 m s −1 , respectively. Therefore, the sounding data actually were more representative for observed dry squalls than for tornadoes.
For the event of 18 June 2017, extremely high convective instability (CAPE ≥ 3000 J kg −1 ) was observed over the eastern part of the Kurgan region between 0900 and 1500 UTC ( Table 2). LLS and DLS reached 11 and 21 m s −1 , respectively. The hodograph ( Figure 3) displayed a pronounced 'kink' near the 900 hPa surface, where a combination of weak directional shear and strong speed shear abruptly changed to strong directional shear, which was likely associated with a low-level mesocyclone [76]. Other parameters were also above their critical levels. SRH exceeded 450 m 2 s −2 , the maximum value of EHI reached 5.0, SCP and STP amounted to 3.5 and 2.9, respectively, while CIN amounted to −72 J kg −1 . Thus, the atmospheric environments were exceptionally favorable for the development of a severe weather outbreak, including strong tornadoes.

Convective Storms Evolution
The evolution of convective storms was analyzed using the Meteosat-8 images with the 15-min time resolution. Weather station data and eye-witness reports were also utilized.
On 3 June 2017, a deep convection started already at 0400 UTC on the south of the Chelyabinsk region. In the next few hours, the MCS, which developed near the occlusion point of the cyclone, increased its size up to 250 km and moved northward; severe wind gusts (up to 27 m s −1 , see Table A3 for details) and large hail (with a diameter up to 35 mm, Table A3) were reported at the weather stations in the Chelyabinsk region.
At 0830 UTC, an isolated convective cell with cloud top temperature (CTT) of −60 • C formed over northeastern Bashkortostan, behind the front. According to eye-witness observations, the cell contained a mesocyclone. At 1115 UTC, a well-pronounced OT with CTT −62 • C was observed by the Meteosat-8 image near Staroutkinsk (the western part of the Sverdlovsk region), where the first IF2 tornado of that day was reported (Figure 4a,b). Next 30 min, OT transformed to the cold-ring feature with a diameter of 80 km (Figure 4c,d). At the same time, the second significant tornado passed near Visim. After 1200 UTC, the storm weakened. Another isolated storm, which probably contained a mesocyclone, formed at 0945 UTC, 100 km to the north of the track of the first storm. It passed between 1100 and 1200 UTC near Biser weather station (WMO ID 28138), where a large hail (32 mm) was reported.
Afterward, between 1130 and 1330 UTC, another convective storm formed on the cold front and moved through the Chelyabinsk region and the western part of the Sverdlovsk region. It caused widespread damaging winds (22-27 m s −1 ) and numerous mostly weak tornadoes (IF0-IF1, one IF2), which were identified by the high-resolution satellite images. It is of note that severe wind gusts caused more substantial damage than tornadoes. After 1330 UTC, the MCS merged with two isolated storms and formed mesoscale convective complex (MCC) with a major axis diameter of about 400 km, and CTT decreased down to −65 • C. The MCC further yielded several weak (IF0-IF1) tornadoes and one IF2-tornado, severe wind gusts, and large hail in the north of the Sverdlovsk region and northeast of the Perm region (Tables A2 and A3). It reached a maximum intensity near Kachkanar, where 400 ha of forest were totally damaged by a windstorm.

The Main Characteristics of Tornadoes and Associated Severe Weather Events
We found 28 tornadoes for the 3 June 2017 event and 9 tornadoes for the 18 June 2017 event ( Figure 6, Table A2) by collecting and analyzing eye-witness observations and damage reports in the media, as well as evaluating forest damage using satellite data. Among all tornadoes, 29 cases were found with high-resolution satellite images only, while 8 tornadoes were reported by eye-witnesses. Information on seven tornadoes was confirmed by two or more independent data sources. We estimated the certainty of events (see [2] for details of this procedure) and defined the high level of certainty for 32 tornadoes, and the medium level for five tornadoes.
Among 28 tornadoes of 3 June 2017, the strongest tornado passed near Visim and damaged forest on the area of 4.97 km 2 (Figure 7b,c). Its intensity substantially depends on the method of estimating (Table A2). In particular, according to the method [66], the intensity is estimated as EF4, since more than 80% of trees were blown down or snapped in 100 m × 100 m forest areas. However, tree debarking (which is considered a signature of ≥IF3 intensity) cannot be confirmed by aerial images; additionally, it is important to note that spruce and birch trees were mostly uprooted than snapped. Therefore, based on the recommendations for the International Fujita Scale [65], the tornado intensity should be estimated at the IF2 level, taking into account the domination of low-firmness trees in the analyzed region. On the contrary, geometrical characteristics of the tornado path, i.e., the length of 33.3 km (the longest tornado path for both outbreaks) and the maximum width of 1590 m (the second widest) put the tornado intensity to the at least F3 category with a 90% probability (Table A2). Considering the three estimates and mismatch between F-/IF-and EF-scales, we speculate that this tornado likely has IF3 intensity.
Another significant tornado that day passed through Staroutkinsk, where it seriously damaged dozens of houses and damaged forest in the area of 1.57 km 2 (Figure 7a). Other tornadoes passed through forested area and caused local windthrows (mostly corresponding to the IF0 or IF1 intensity levels). The total area of stand-replacing windthrow caused by the severe weather outbreak of 3 June 2017 was estimated as 15.4 km 2 .

The Main Characteristics of Tornadoes and Associated Severe Weather Events
We found 28 tornadoes for the 3 June 2017 event and 9 tornadoes for the 18 June 2017 event ( Figure 6, Table A2) by collecting and analyzing eye-witness observations and damage reports in the media, as well as evaluating forest damage using satellite data. Among all tornadoes, 29 cases were found with high-resolution satellite images only, while 8 tornadoes were reported by eye-witnesses. Information on seven tornadoes was confirmed by two or more independent data sources. We estimated the certainty of events (see [2] for details of this procedure) and defined the high level of certainty for 32 tornadoes, and the medium level for five tornadoes.
Among 28 tornadoes of 3 June 2017, the strongest tornado passed near Visim and damaged forest on the area of 4.97 km 2 (Figure 7b,c). Its intensity substantially depends on the method of estimating (Table A2). In particular, according to the method [66], the intensity is estimated as EF4, since more than 80% of trees were blown down or snapped in 100 m × 100 m forest areas. However, tree debarking (which is considered a signature of ≥IF3 intensity) cannot be confirmed by aerial images; additionally, it is important to note that spruce and birch trees were mostly uprooted than snapped. Therefore, based on the recommendations for the International Fujita Scale [65], the tornado intensity should be estimated at the IF2 level, taking into account the domination of low-firmness trees in the analyzed region. On the contrary, geometrical characteristics of the tornado path, i.e., the length of 33.3 km (the longest tornado path for both outbreaks) and the maximum width of 1590 m (the second widest) put the tornado intensity to the at least F3 category with a 90% probability (Table A2). Considering the three estimates and mismatch between F-/IF-and EF-scales, we speculate that this tornado likely has IF3 intensity. Most of the tornadoes of 18 June 2017 passed mainly over a treeless area ( Figure 6) and were reported by eye-witnesses. The strongest tornado hit Maloye Pesyanovo village at 1135 UTC and blew down four strong log-houses; also, tree debarking was reported (Figure 8c). According to the recommendations for the International Fujita scale [65], the near complete destruction of log houses can be considered as an indicator of IF4 damage since the sturdiness of such houses is somewhat equal to sturdiness of 'strong framehouses' (Thilo Kühne, personal communication). The violent character of this tornado is also supported by two other intensity estimates. Particularly, the tornado has the EF5 intensity level according to the method [66] (up to 100% of trees were blown down or snapped in 100 m × 100 m forest areas). Given tornado path length (20.3 km) and maximum width (1750 m, the widest for both outbreaks), the minimal tornado intensity was estimated as F3 with 90% probability. Consequently, our estimates confirm the formation of the first IF4 tornado in Russia in the 21st century (after the previous one in 1984 [10]) and the first-ever IF4 tornado reported beyond the Urals. Another significant tornado that day passed through Staroutkinsk, where it seriously damaged dozens of houses and damaged forest in the area of 1.57 km 2 (Figure 7a). Other tornadoes passed through forested area and caused local windthrows (mostly corresponding to the IF0 or IF1 intensity levels). The total area of stand-replacing windthrow caused by the severe weather outbreak of 3 June 2017 was estimated as 15.4 km 2 .
Most of the tornadoes of 18 June 2017 passed mainly over a treeless area ( Figure 6) and were reported by eye-witnesses. The strongest tornado hit Maloye Pesyanovo village at 1135 UTC and blew down four strong log-houses; also, tree debarking was reported (Figure 8c). According to the recommendations for the International Fujita scale [65], the near complete destruction of log houses can be considered as an indicator of IF4 damage since the sturdiness of such houses is somewhat equal to sturdiness of 'strong framehouses' (Thilo Kühne, personal communication). The violent character of this tornado is also supported by two other intensity estimates. Particularly, the tornado has the EF5 intensity level according to the method [66] (up to 100% of trees were blown down or snapped in 100 m × 100 m forest areas). Given tornado path length (20.3 km) and maximum width (1750 m, the widest for both outbreaks), the minimal tornado intensity was estimated as F3 with 90% probability. Consequently, our estimates confirm the formation of the first IF4 tornado in Russia in the 21st century (after the previous one in 1984 [10]) and the first-ever IF4 tornado reported beyond the Urals. Atmosphere 2020, 11, x FOR PEER REVIEW 15 of 35

Results of Modeling
In this section, we present the results of cloud-resolving simulation of tornadic thunderstorms using the WRF-ARW mesoscale model [39]. To estimate the forecast accuracy, we compared the simulated severe storms, i.e., their paths and intensity as well as cloud top characteristics, with the observed tornado tracks and the Meteosat-8 images. In Section 4.1, we describe the results of the performed experiments with various model grid spacings and with the 12-h lead time. Section 4.2 contains an assessment of the model performance with one selected grid spacing but with various (24-and 36-h) lead times.

Determination of the Most Appropriate Model Grid Spacing
We performed three simulations for both tornado events, with model settings described in Section 2.3 with a 12-h lead time (model start from 0000 UTC 3 June 2017 and 0000 UTC 18 June 2017). To assess the forecast reliability, firstly we determined the presence of potentially tornadic storms in the model output (see Section 2.3 for details). Then, if a storm was predicted by the model, we estimated the location and timing errors by comparing the simulated storm location with that observed on the Meteosat-8 images, and simulated the tornadic storm track with the Landsat-based tornado tracks. Figures 9 and 10 and Table 3

Results of Modeling
In this section, we present the results of cloud-resolving simulation of tornadic thunderstorms using the WRF-ARW mesoscale model [39]. To estimate the forecast accuracy, we compared the simulated severe storms, i.e., their paths and intensity as well as cloud top characteristics, with the observed tornado tracks and the Meteosat-8 images. In Section 4.1, we describe the results of the performed experiments with various model grid spacings and with the 12-h lead time. Section 4.2 contains an assessment of the model performance with one selected grid spacing but with various (24-and 36-h) lead times.

Determination of the Most Appropriate Model Grid Spacing
We performed three simulations for both tornado events, with model settings described in Section 2.3 with a 12-h lead time (model start from 0000 UTC 3 June 2017 and 0000 UTC 18 June 2017). To assess the forecast reliability, firstly we determined the presence of potentially tornadic storms in the model output (see Section 2.3 for details). Then, if a storm was predicted by the model, we estimated the location and timing errors by comparing the simulated storm location with that observed on the Meteosat-8 images, and simulated the tornadic storm track with the Landsat-based tornado tracks. Figures 9 and 10 and Table 3      The tornadic storm of 3 June 2017 was reproduced by the model only with the version with the 3-km grid spacing and no nested grid (Figure 9a). An isolated storm (a meso-scale low) with SRH up to 1075 m 2 s −2 and composite reflectivity~58 dBz was predicted between 1200 and 1400 UTC, reaching the maximum intensity at 1300 UTC according to the model output. However, the simulated intensity of the severe storm was underestimated since wind gusts did not exceed 17 m s −1 . The timing error of the forecast was about 1-1.5 h.

WRF-Simulated Storm Parameters (Maximum Values in 50-km Radius around a Tornado
The WRF model configuration with a 7.2-km grid resolution did not predict isolated storms that could generate tornadoes (Figure 9b). The simulation with a 9-km grid resolution and a 3-km nested grid (Figure 9c) predicted two storms with high composite reflectivity (>50 dBz) that moved east of the observed tornado track and had no above-mentioned storm signatures. We failed to simulate strong wind gusts (>20 m s −1 ) with any WRF model settings. In general, we can conclude that none of the three simulations performed at 00 UTC 3 June 2017 allowed prediction of the tornadoes that occurred on 3 June 2017.
The 18 June 2017 tornadic storm over the Kurgan region was simulated by the WRF model accurately, in particular with the version with the 3-km grid spacing and no nested grid (Figure 10a). At 1300 UTC, the model predicted well-pronounced meso-scale low severe convective storm with composite reflectivity up to 64 dBz, severe wind gust up to 31 m s −1 , and SRH~1200 m 2 s −2 . The simulated severe storm track and observed tornado track had a displacement less than 10 km. Considering that the model with the 3-km resolution does not allow simulation of the tornado itself, this forecast can be estimated as successful.
However, the simulated severe storm passed near Maloye Pesyanovo village around 1300 UTC, whereas the actual tornado occurred at 1135 UTC (Table A2). This time lag (1 h 25 min) is presumably associated with the GFS model assimilated data, since the lag time is independent of the WRF model configuration.
The forecasts of tornadic storm of 18 June 2017 with other WRF model settings (Figure 10 b,c) were less accurate. The model predicted a severe storm track 30 km to the east of the observed tornado track. However, this storm was not confirmed neither by the Meteosat-8 images nor the forest damage analysis. Instead of a tornadic storm, the model predicted a squall line, which in addition was slightly shifted to the west compared to the tornado track.
For both cases, the 10-30-km displacement between the simulated storm tracks and observed tornado tracks is generally in line with results obtained previously [22,42,77]. Our findings on the simulation of two or more tornadic storms instead of one observed storm (as in the case of 18 June 2017) is in agreement with [22,45].
Comparing three WRF model configurations for both events, we found the most appropriate accuracy of tornadic storm forecast using the WRF model configuration with the 3-km grid spacing and no nested grid. Further, we considered this model configuration as optimal for a short-term forecast of severe storms for the analyzed events.

Validation of the Simulation Results with the Meteosat-8 Images
The MCSs characteristics (CTT and composite reflectivity) simulated with all three WRF model configurations were overlapped with the Meteosat-8 images. Table 4 presents the simulated and satellite-derived CTTs, as well as location and timing errors that were estimated by comparing simulated and observed spatial position of convective storms. Figures 11 and 12 present the convective storms at the time of tornado occurrence based on the Meteosat-8 images and the WRF simulations with the optimal model configuration (3-km resolution, no nested grid).

Sensitivity of Simulation Results to Forecast Lead Time
The WRF model with a 3-km resolution and no nested grid started at 0000/1200 UTC 2 June 2017 and 17 June 2017 and was used to estimate the influence of the forecast lead time on the accuracy of simulation. Table 5 summarizes the obtained results, and Table 6 presents the comparison of simulated MCSs characteristics with the Meteosat-8 data.
The model predicted the convective storms with high composite reflectivity and strong wind gusts independently of the forecast lead time. The most reliable forecast of the event of 3 June 2017 was obtained from 1200 UTC 2 June 2017 (i.e., with the 24-h lead time) when the model predicted a storm with a well-pronounced meso-scale low-pressure system, very high SRH (up to 1350 m 2 s −2 ), and wind gusts up to 28 m s −1 (Figure 13b). The forecast with the 36-h lead time is also relatively successful, since the model predicted the convective storm but shifted 50 km to the east (Figure 13b). Both forecasts (24 and 36 h) have a substantial timing error (1.25-1.5 h). In particular, at 1200 UTC, the simulated storm was 20 km to the south of the start of the first tornado track. In fact, at this time, the tornado had already passed near Visim town, 70 km to the north of its simulated location.
In the second case (18 June 2017), the model predicted the convective storm with strong wind gusts (up to 30 m s −1 ) but with relatively low values of SRH (≤600 m 2 s −2 ). Thus, the forecast of the 18 For the first case (3 June 2017, Figure 11), one can see that the WRF-simulated and satellite-observed location and geometry of the northern severe storm match rather well. Simulated and observed CTT values range from −55 to −60 • C. It is of note that the cold-ring signature on the cloud top is also well-reproduced by the model. Nevertheless, the forecast cannot be considered as successful, since the model entirely missed the main MCS that formed on the cold front 100 km to the south of the first storm and caused wind-related damage.
For the second case (18 June 2017, Figure 12), the tornado-generating MCC with a diameter exceeding 400 km was well-predicted by the model. Several well-detected OTs and cold-ring signatures are seen on the cloud top. It is important to note that the OTs correspond to the zones with simulated high composite reflectivity. One OT moved along the tornado track with only a 10-km displacement. The main shortcoming of the simulation is the substantial time error. Indeed, despite the simulated and observed MCC tracks having similarity, the simulated MCC moved with a 1.5-h lag compare to the satellite-based one (Table 4).

Sensitivity of Simulation Results to Forecast Lead Time
The WRF model with a 3-km resolution and no nested grid started at 0000/1200 UTC 2 June 2017 and 17 June 2017 and was used to estimate the influence of the forecast lead time on the accuracy of simulation. Table 5 summarizes the obtained results, and Table 6 presents the comparison of simulated MCSs characteristics with the Meteosat-8 data.  The model predicted the convective storms with high composite reflectivity and strong wind gusts independently of the forecast lead time. The most reliable forecast of the event of 3 June 2017 was obtained from 1200 UTC 2 June 2017 (i.e., with the 24-h lead time) when the model predicted a storm with a well-pronounced meso-scale low-pressure system, very high SRH (up to 1350 m 2 s −2 ), and wind gusts up to 28 m s −1 (Figure 13b). The forecast with the 36-h lead time is also relatively successful, since the model predicted the convective storm but shifted 50 km to the east (Figure 13b). Both forecasts (24 and 36 h) have a substantial timing error (1.25-1.5 h). In particular, at 1200 UTC, the simulated storm was 20 km to the south of the start of the first tornado track. In fact, at this time, the tornado had already passed near Visim town, 70 km to the north of its simulated location.
In the second case (18 June 2017), the model predicted the convective storm with strong wind gusts (up to 30 m s −1 ) but with relatively low values of SRH (≤600 m 2 s −2 ). Thus, the forecast of the 18 June 2017 storm with the 24-and 36-h lead time in principal can be used to warn severe wind gusts but not tornadoes. The forecast time error was approximately the same as for the 12-h lead time forecast.
Therefore, we showed that the WRF model forecasts obtained with 24-and 36-h lead times can also be appropriate to predict the threat of severe storms and tornadoes during the analyzed events. This is in line with previously published results for the U.S. territory [29,33]. Atmosphere 2020, 11, x FOR PEER REVIEW 24 of 35

Discussion and Conclusion
In this study, we evaluated the environments of the formation of two tornado outbreaks that occurred east of the Ural Mountains on 3 June and 18 June 2017 using different data. For tornado diagnostics, we used news reports, eyewitness reports, ground-based meteorological observations, sounding data, global NWP models data, synoptic charts, satellite data on forest disturbances, and data of the specially conducted aerial imaging. We also estimated the accuracy of short-term forecasting of outbreaks with the WRF-ARW atmospheric model.
Both analyzed outbreaks caused substantial damage associated with tornadoes, severe wind gusts, and large hail. We found the 3 June 2017 outbreak is rather rare in terms of the number of tornadoes-we confirmed the formation of 28 tornadoes. Four of them were significant (≥IF2 intensity), including one strong tornado (IF3). Convective storms moved mainly over forest-covered area, which made it possible to identify tornado damage tracks in forests using the Sentinel-2 data, high-resolution satellite images, and data of the aerial imaging performed for the strongest tornado

Discussion and Conclusions
In this study, we evaluated the environments of the formation of two tornado outbreaks that occurred east of the Ural Mountains on 3 June and 18 June 2017 using different data. For tornado diagnostics, we used news reports, eyewitness reports, ground-based meteorological observations, sounding data, global NWP models data, synoptic charts, satellite data on forest disturbances, and data of the specially conducted aerial imaging. We also estimated the accuracy of short-term forecasting of outbreaks with the WRF-ARW atmospheric model.
Both analyzed outbreaks caused substantial damage associated with tornadoes, severe wind gusts, and large hail. We found the 3 June 2017 outbreak is rather rare in terms of the number of tornadoes-we confirmed the formation of 28 tornadoes. Four of them were significant (≥IF2 intensity), including one strong tornado (IF3). Convective storms moved mainly over forest-covered area, which made it possible to identify tornado damage tracks in forests using the Sentinel-2 data, high-resolution satellite images, and data of the aerial imaging performed for the strongest tornado of this outbreak. The characteristics of tornado paths were determined, i.e., tornado width and length. On the contrary, the outbreak 18 June 2017 occurred over mainly treeless area. Therefore, it is probably that the number of tornadoes within this outbreak is underestimated, since they cannot be found with satellite data. In addition, the geometric characteristics of tornado paths were determined not for all tornadoes. The outbreak includes four significant tornadoes, including one violent tornado that hit Maloye Pesyanovo village-the first IF4 tornado in Russia in the 21st century and the first-ever IF4 tornado recorded beyond the Ural Mountains.
Only four tornadoes that passed through settlements and one tornado that was covered by aerial imaging have solid estimates of tornado intensity. The intensity of other tornadoes was determined rather approximately. We used three different methods to estimate their intensity based on satellite data only. We used the recommendations for International Fujita scale [65], which allowed us to discriminate tornadoes with IF0 intensity among all others. Using satellite data, we failed to reveal the signatures of trees' damage indicating higher tornado intensity; particularly, it was impossible to see whether trees were snapped or debarked. We also used our method developed in [10,62], which is based on the relationship between Fujita-scale tornado intensity and tornado path width and length [68]. This method provides only a minimal intensity rating (i.e., ≥F2) and cannot be used when there is a lack of data on tornado track characteristics. In addition, the tornado intensity estimate may be incorrect for short-lived tornadoes. The third used approach is based on the method of Godfrey and Peterson [66], who tied EF-scale tornado intensity with forest damage severity. We found this method generally overestimated the intensity of the analyzed tornadoes. It is of note that the method had been developed using linkages between forest wind resistance and wind speed obtained for North American forest species [66]. Presumably, this linkage is not well suited for northern Eurasian forests. This issue deserves further research.
The analysis of synoptic-scale environments revealed the similarity of the two outbreaks. Both events formed in the afternoon hours (at local time), when the development of a deep convection reached a daily maximum [2]. Additionally, both events formed near the cold polar front in the warm part of deepening southern cyclones, close to their occlusion points. Tropical air masses spread from Kazakhstan, and the high precipitable water content (up to 42 mm on 18 June 2017 and 30 mm on 3 June 2017) was related mainly to strong convergence near the cyclone center and local evaporation from the well-moistened underlying surface. Such synoptic-scale environments are favorable for tornado occurrence both in the Ural region [11,16] and in the European part of Russia [6,14,16,19]. The synoptic situations of both events correspond to the second typical scheme of tornado formation over northern Eurasia according to Snitkovsky [4]. Based on the statistics for the European part of Russia, Snitkovsky had found that this type yields mostly weak tornadoes [4]. However, events analyzed in the study resulted in IF3 and IF4 tornadoes, which indicates the need to reconsider Snitkvosky schemes for other Russian regions (i.e., the Ural, Siberia, Far East). This can be achieved using data from the new climatology presented in [2].
The characteristics of tornado-generating MCSs were substantially different in the first and second cases. On 3 June 2017, two significant tornadoes were induced by the long-lived severe storm that was formed in the warm sector. Additionally, numerous weak tornadoes occurred when there was a mesoscale convective system on the cold front. On 18 June 2017, tornadoes were associated with mesocyclones embedded into a rapidly increasing mesoscale convective complex that formed on the cold front. In both cases, a well-pronounced OT and cold-ring signatures on the cloud top coincided in space and time with tornado occurrence.
Short-term forecast based on the global NWP models (GFS and GEM) produced critical values of diagnostic variables (severe-weather indices) that indicate tornado formation for both events. In particular, a rare combination of strong convective instability (CAPE > 1500 J/kg) and low-level wind shear (>10 m s −1 ) was found. Such environments are favorable for tornado occurrence in the U.S. [21,51] and in Europe [17,59]. An exceptionally high value of EHI (up to 4.0 and 5.0 for the first and second events) and STP (2.5 and 2.9) also unambiguously pointed to a high risk of formation of strong tornadoes.
The specially conducted short-term forecasts with the convection-permitting WRF model were rather successful in both cases. The model reproduced the formation of severe convective storms, which somewhat coincides in space and time with the observed tornado tracks. However, for the 3 June event, the model failed to reproduce the main storm (the second storm of that day), which induced widespread severe wind. For the 18 June event, the model reproduced the formation of the MCC with several well-pronounced meso-scale lows and wind gusts > 30 m s −1 . In both cases, the WRF model configuration with the 3-km grid spacing and no nested grid provided the most reliable forecast. We also evaluated the forecast sensitivity to the lead time (12-36 h) and found that the model simulated the convective storms somewhat independently of the forecast lead time. We should stress that it is necessary to evaluate more cases of tornadoes and tornado outbreaks to determine the optimal configuration of the WRF model and estimate the sensitivity of simulated storm tracks characteristics (e.g., their location and timing errors) to a forecast lead time.
Operational short-term tornado forecasting is not implemented in Russia. Our analysis demonstrates that combined use of the global NWP model outputs indicating high values of severe-weather indices and the WRF model forecast outputs explicitly indicating convective storms and meso-cyclones formation could be used to predict the high probability of strong tornado occurrence. For both analyzed events, such tornado warning forecast could help local authorities to take early action on population protection. However, it should be highlighted that the frequency of false alarms in WRF-based short-term forecasts of tornadic storms remains unexplored, at least for the territory of northern Eurasia. To obtain such estimates, one needs to build and analyze long-term series of mesoscale NWP models output data, e.g., The Spring Experiment conducted in the USA [24,28]. In the future, we are planning to set-up and analyze such experiments for several Russian regions.

Conflicts of Interest:
The authors declare no conflict of interest.   [3] The following notations are used: LFC-height of level of free convection, EL-height of equilibrium level, T v,parcel -virtual temperature of the individual parcel (lifted from surface), T v,env -virtual temperature of the environment, T' 500 -virtual temperature at 500 hPa, T' p,500 -virtual temperature of lifted parcel at 500 hPa, V-wind speed, k-vertical unit vector, C-moving velocity of cell, V h -wind vector, MUCAPE-most unstable CAPE, SBLCL-Surface-based Lifted Condensation Level. Subscripts denote vertical levels (surf-near surface; 900 and 500-corresponding isobaric levels; 0-3-the layer from surface to 3 km height; 0-1-the layer from the surface to 1 km height).    [65], the second value obtained based on forest damage analysis according to [64], the third value obtained based on tornado path length and width according to [10,62].