Observations and modelling of the winter thunderstorm on 4 February 2022 at the Milešovka meteorological observatory

The study analyses a winter thunderstorm that passed over the Milešovka meteorological observatory on 4 February 2022, between 2300 and 2330 UTC. Lightning was recorded directly over the observatory by both the observer and the EUCLID lightning network at 2320 UTC. To analyse the state of the atmosphere at the time when the lightning occurred, we used data from the X‐band Doppler polarimetric radar and the Ka‐band Doppler polarimetric vertical profiler, both located at the observatory. We also applied data from the Meteosat Second Generation satellite, and data from standard meteorological instruments located at the observatory. In addition, we run our cloud electrification model to simulate cloud electrification of the winter thunderstorm to find out whether the model develops conditions suitable for the occurrence of lightning and if so, under what circumstances. Our results show that the lightning appeared at the very end of the storm passage defined by high radar reflectivity. At the same time, it is clear from the radar observations that before lightning occurred, the cloud contained hydrometeors (graupel, cloud or rain water, and ice or snow) which are commonly associated with charge separation by collisions. Our analysis of the radar data also suggests that in at least several parts of the cloud the electric field was strong. Although the cloud top height was very low compared to summer storms, the model results indicate conditions suitable for lightning occurrence. However, uncertainty remains on how to properly formulate the initial conditions for model simulations for this type of storm which was shallow and occurs rarely in winter.

studies on winter storms can be found in literature.Most of the latter focused on the statistical features such as spatial distribution and frequency of winter storms.For instance, March et al. (2016) presented winter lightning maps of selected regions in the Northern Hemisphere.They defined four different degrees of winter lightning activity and prepared a guidance on risk assessment for tall structures and wind turbines.Kolmašová et al. (2022) used the World Wide Lightning Location Network (WWLLN; http://wwlln.net/)data to investigate lightning strokes that hit Northern Europe during the unusually stormy winter of 2014/2015.Morgenstern et al. (2022) examined and compared meteorological processes responsible for lightning in winter with those responsible for lightning in summer.They used cluster analysis and principal component analysis to find physically meaningful groups of processes in ERA5 atmospheric reanalysis data and lightning data for northern Germany.They found two types of thunderstorms: (i) wind-field thunderstorms and (ii) convective available potential energy (CAPE) thunderstorms.
Winter thunderstorms were also investigated in terms of the connection of lightning occurrence with other phenomena.For example, Pizzuti et al. (2022) studied lightning occurrence including the potential of superbolts (i.e.short-lived optical phenomena above thunderstorms) and possible causes in the area of the English Channel.Several studies focused on lightning strokes associated with transient luminous events, for example, Takahashi et al. (2003), Hayakawa et al. (2004), Yair et al. (2015), and Peterson and Kirkland (2020).Finally, there are published studies which analysed the structure of winter storms and investigated their evolution, for example, Altaratz et al. (2021).Altaratz et al. (2021) found out that the lightning activity corresponds to the mature stage of the thunderstorm with the maximum frequency of flashes at the end of this stage because in mature stage, cells are developed upwards in the updraught which makes the radar reflectivity intensify and push it to higher elevations.
In contrast, very little work was focused on modelling the evolution of the electric field in winter thunderstorms as compared to that in summer thunderstorms.Fierro et al. (2013) studied a continental winter storm using an electrification model based on bulk lightning approach with explicit charging and discharging physics.They implemented it into the Weather Research and Forecasting (WRF) numerical weather prediction (NWP) model.They ran the model with a horizontal mesh size of 3 km and the model did not generate any lightning during the forecast, which was consistent with observations.The electric field was found several orders of magnitude weaker in the case of the winter storm than that of typical summer storms.Another study which referred to a simple one-dimensional (1D) mathematical model to simulate evolution of a winter storm was that by Altaratz et al. (2021) mentioned above.
Winter thunderstorms are very rare in Central Europe and if they occur they are usually related to passages of pronounced cold fronts (Munzar and Franc, 2003), which is the case of our study.This study focuses on a winter thunderstorm which occurred on 4 February 2022, and hit the Milešovka meteorological observatory, which is located in Central Europe.We analyse the structure of the storm using both conventional meteorological data as well as radar data and satellite data.We take the advantage of two radars for this study: (i) X-band Doppler polarimetric radar, and (ii) Ka-band Doppler polarimetric vertical profiler, both located at the observatory.In addition, we use an idealized cloud electrification model to simulate the winter thunderstorm.We investigate whether and under what circumstances the model develops conditions suitable for lightning.To initialize the model, we use prognostic data from the ICON-EU (https://code.mpimet.mpg.de/projects/iconpublic) NWP model.
The study has two main objectives.The first objective is to determine if the observed data, especially radar data, have specific features that can be related to cloud electrification and discharges.The second objective is to determine under what conditions a simplified cloud model containing a description of cloud electrification is able to develop an electric field that leads to discharges in wintertime conditions.These two objectives of the study are related, but it should be noted that the second objective does not aim at predicting and comparing the model evolution of the cold front including lightning activity with observations.The cloud model with electrification that we use cannot predict accurately the development of a cold front and certainly cannot accurately predict the occurrence of lightning in time and space.The modelling part of the study answers the question of whether and under what conditions a winter thunderstorm with lightning develops based on data derived from real measurements.Partly, it also aims at verifying how the process of cloud electrification and discharge generation depend on the horizontal and vertical resolution of the model.Increasing the model resolution will affect the model fields, but it is questionable how this will affect the discharge generation, because discharge generation is a strongly nonlinear process that depends on overpassing thresholds of various variables.
The article is organized as follows.Localization and equipment of the Milešovka meteorological observatory is described in Section 2. Section 3 deals with lightning data, while Section 4 presents basic meteorological data.Sections 5 and 6 show our results.Section 5 describes the vertical structure of the winter thunderstorm derived from data of the two radars, whereas Section 6 focuses on the application of our cloud electrification model to simulate

MILEŠOVKA METEOROLOGICAL OBSERVATORY
The Milešovka meteorological observatory (50 • 33 ′ 17 ′′ N, 13 • 55 ′ 57 ′′ E, 837 m a.s.l., WMO No. 11464) is located in the northern part of the Czech Republic at the top of the Milešovka mountain (Figure 1).Milešovka is a pointed hill exceeding its surroundings by more than 300 m vertically.Meteorological measurements at the observatory are thus close to those of the free atmosphere above the boundary layer.About 20 km away to the north and north-northwest, there is a mountain ridge elevated at 1,000 m approximately, which often influences the air characteristics, for example, wind speed, temperature, and humidity at the Milešovka observatory.
The observatory is part of the national and international meteorological networks.It provides routine and specific meteorological and climatological measurements with a 24/7 service.The station's 360 • view as well as its long-term series of measurements (since 1905;Štekl and Podzimek, 1993) make the observatory unique in the European context.
The observatory is equipped with standard meteorological and climatological instruments (measurements of pressure, air temperature, soil temperature, humidity, wind direction and wind speed, global radiation, and precipitation by several devices).The instruments also include a Vaisala ceilometer CL51 and Thies Laser Precipitation Monitor disdrometer.Further, the Milešovka observatory is equipped with two research meteorological radars: (i) a vertically pointing cloud radar MIRA 35c, and (ii) an X-band weather radar FURUNO.Their description follows below.

Cloud radar MIRA 35c
The vertically pointing cloud radar MIRA 35c working at 35 GHz (Ka-band) was installed at the Milešovka observatory in 2018.It is a Doppler polarimetric radar operating with a vertical range of 28.9 m and a time step of 2 s.It can measure up to 14 km vertically.More radar parameters can be found in Sokol et al. (2018).
The basic quantities measured by the radar are: Doppler velocity, power in co-and cross-channels, and phase in co-and cross-channels of reflected signal.To analyse structure of the winter thunderstorm, we use equivalent radar reflectivity (Ze), Doppler velocity (V), Linear Depolarization Ratio (LDR), and Co-polar correlation coefficient (RHO).All these quantities are derived from the basic measured quantities.We also estimate the vertical air velocity (Va) and we classify hydrometeors into six types (rain, hail, water droplets, graupel, snow and ice) and their combinations, based on radar data.In this study, we use the HClass classification algorithm to distinguish hydrometeors, as proposed by Sokol et al. (2018) and refined in Sokol et al. (2020).Calculation of Va requires wind de-aliasing.We applied the de-aliasing procedure presented by Sokol et al. (2020) to estimate Va.

X-band weather radar FURUNO WR2120
The X-band weather radar FURUNO WR2120 is a Doppler fully polarimetric X-band weather radar which was installed at the Milešovka observatory in 2020.Basic operational characteristics of the radar are the following: horizontal resolution of 150 m, frequency of 9.4 GHz, antenna beam-width of 2.5 • , and horizontal range of 35 km.More details about the radar and its configuration can be found in Bobotová et al. (2022).
More details about the standard radar variables listed above can be found, in, for example, Ryzhkov and Zrni ć (2019).We use the measured data to classify hydrometeors.We classify seven types of hydrometeors by applying the XClass algorithm on the data, which was developed for the radar by Bobotová et al. (2022).
To investigate the cloud structure above the Milešovka observatory, we apply a scanning strategy emphasizing vertical cross-sections.The radar scanning strategy consists of two steps: (i) the radar performs 7 horizontal PPI (plan position indicator) clockwise scans for elevation angles 1.1 • , 1.7 • , 2.5 • , 4 • , 6 • , 10 • and 25 • ; (ii) the radar makes 6 RHI (range height indicator) scans, that is, 6 vertical cross-section scans for elevation angles from 3 • to 90 • .The difference between adjacent elevation angles is non-equidistant but not larger than 0.5 • .The azimuths of the scans are set in a way that they can be folded into three full cross-sections.The first is oriented from south to north, the second is rotated 240 • clockwise to the first, and the third is rotated 300 • to the first.One cycle of all the scans lasts 170 s.

LIGHTNING DATA
The winter thunderstorm, which was recorded at the Milešovka meteorological observatory on 4 February 2022, was accompanied by lightning flashes detected at 23:20:08 UTC by several sources.We have an eyewitness testimony from an observer who had a position at the Milešovka station and wrote down in the book of records that there was an electric discharge in the immediate vicinity of the station.Two cloud-to-cloud discharges were recorded by the EUCLID (European Cooperation for Lightning Detection, https://www.euclid.org/)network at that time at a distance of less than 200 m from the observatory.In addition, the EUCLID also recorded several lightning discharges 3 min later at a distance of 3-8 km in the southeast direction from the station.Several discharges were also recorded close to the observatory at about the same time by the LINET lightning network (https://www.nowcast.de/en/solutions/linet-systems/) and the Blitzortung network (https://www.blitzortung.org).Some studies have shown that electrical discharges can be triggered by tall structures such as skyscrapers or telecommunication towers, or by gravity waves that are affected by rugged orography (Bech et al., 2013;Warner et al., 2014).From this viewpoint the shape of the hill and the shape of the observatory including a telecommunication tower located near the observatory could have a triggering effect, but we cannot confirm it for the winter storm investigated in this article.The hilly terrain around the Milešovka hill can also be a source of gravity waves; however, given that lightning was observed on the cold front in a flat area even before the front entered the mountainous area, we believe that the orography did not play an essential role in lightning generation.

METEOROLOGICAL CONDITIONS OF THE WINTER THUNDERSTORM ON FEBRUARY 2022
A cold front crossed the northern part of the Czech Republic during the night from 4 February to 5 February 2022.There was a strong northwesterly airflow resulting from a low situated above the Norwegian Sea and a high situated north of the Azores.The front was characterized by a narrow band of clouds that crossed the Milešovka observatory between 2300 and 2330 UTC approximately.As the cold front passed, individual lightning flashes appeared.The direction of the front movement can be seen in Figure 2a,b which show the radar reflectivity Zh measured by the FURUNO radar at a CAPPI (constant altitude plan position indicator) level of 2 km ASL at 23:17:26 and 23:20:16 UTC (the time of the beginning of the PPI scanning cycle) corresponding mostly to the time of the recorded flashes.The Figure 2c,d show Rhohv and Zdr at CAPPI 2 km ASL.Radar data show that the field of high reflectivity values was not uniform and the measured reflectivity changed rapidly in time and space.It is worth noting that there is a wedge in the east direction showing different values in the Rhohv and Zdr values which is not very clear in the Zh values.This wedge is also absent in other variables (not depicted), so we assume that it should not be any artefact caused by temporal shadowing or technical problems of the radar.
The passage of the cold front is evident in the evolution of surface values measured at the observatory (Figure 3).Depicted values of surface pressure, wind speed and direction, temperature, humidity and cloud base clearly illustrate changes of characteristics of the air masses from 2300 UTC to 0000 UTC.The temperature drop was almost 2 • C within 20 min and the relative humidity increased to 99%.The station was covered by clouds at the time of the discharge (cloud base height indicated by the ceilometer was about 30 m above the station only).The wind speed decreased initially and then slightly increased (after 2320 UTC), and at the same time wind direction changed by 70 • within 20 min.At the end of the cloud belt passage, that is, at the time of the observed discharge, the atmospheric pressure increased by 1 hPa.The measured precipitation was relatively low.Precipitation started around 2300 UTC and reached its maximum after the front passage at 2330 UTC (i.e. 10 min after the lightning discharge).The disdrometer recorded the maximum of precipitation at the same time and it recognized rain and graupel.
Figure 4 shows the data observed by the Meteosat Second Generation (MSG) geostationary satellite.Specifically, it depicts the brightness temperature in the infrared (IR) channel 10.8 μm at 2320 UTC.It confirms the existence of a cloud belt exceeding in height the clouds of its surroundings.The brightness temperature reached about 250 K above the Milešovka observatory, which corresponds to a height of 3,000 m above the surface.When compared to the temperature sounding measurements at the Praha-Libuš station (80 km from the Milešovka observatory) at 0000 UTC on 5 February 2022, the cloud top height corresponded to the height of temperature inversion.This inversion limited the cloud top height to be low.

VERTICAL STRUCTURE OF THE THUNDERSTORM DERIVED FROM MIRA 35C AND FURUNO WR2120 RADAR DATA
Vertical cross-sections measured by FURUNO WR2120 (Figures 5-7) and time sequence of vertical profiles by MIRA 35c (Figures 8 and 9) radars were used to study the vertical structure of the winter thunderstorm.The directions of the cross-sections were predefined.They correspond to angles 0 • , 240 • and 300 • .Note that comparison of the results given by MIRA 35c with those given by FURUNO WR2120 is difficult because the quantities measured and the way they are measured are different for the two radars.MIRA 35c gives data in the vertical direction above the radar only, but with significantly higher temporal and spatial resolutions.It provides more detailed data because the entire Fourier spectrum of the received signal is analysed.On the other hand, FURUNO WR2120 provides the "spatial" view on the horizontal structure and extent of the storm.
Figure 5 shows the RHI cross-sections of measured Zh by the X-band radar at a distance of up to 5 km from the radar located at the point (0,0).The cross-section directions as well as approximate direction of the front movement are displayed in the upper-right panel of Figure 5.The time of the depicted measurement 23:19:16 UTC corresponds to the beginning of the RHI scans and is close to the observed lightning at 23:20:08 UTC.
Figure 6 shows the RHI NW-SE (300 • ) cross-sections of Zh, Zdr, V, W, Rhohv and Kdp quantities measured by the X-band radar.This direction is chosen to best match the direction of the storm motion.Zh values correspond to the lowest Zh values of the wedge in the eastern area of the radar in Figure 2. The wedge direction roughly corresponds to the direction of the cross-sections in Figure 6.The radial velocity (V) in Figure 6 confirms that the flow direction was from the northwest to the southeast, which agrees well with the cold front motion.However, there are clear inhomogeneities in values of V at heights of 2,500 m and above.Sign changes indicate rotational wind motions.
Similar inhomogeneities in values are evident in W and Zdr fields (Figure 6).Larger values of W indicate the presence of various hydrometeors with different terminal velocities, especially that of graupel.Inhomogeneities in values of V, W and Zdr fields at 2,500 m and above correspond to the region where Zh was less than 0 dBZ. Figure 6b shows that there is a strong contrast between positive and negative values in Zdr above 2,550 m in the range from −1 to 1 km from the radar position approximately.These inhomogeneities in values are not consistent with those visible in V and W variables.Some of them are in the region where Zh is from 0 to 5 dBZ, but the majority of them can be found in areas where Zh is very low.The explanation of the variability in Zdr values is not straightforward, because, among other things, the measured values of Zdr and Zh depend on the shape of the hydrometeors, on their orientation, and on the angle at which they are measured by the radar.Since the Zdr field structure is predominantly vertical, we suppose that the direction of the radar measurements does not significantly affect the Zdr values in the range from −1 to 1 km from the radar position approximately.
The inhomogeneities in values of Zdr can indicate that the orientation of ice particles, for example snowflakes, with irregular shapes is disturbed by the airflow in these regions.Positive Zdr indicates horizontally oriented flakes, while negative Zdr indicates vertically oriented flakes.The orientation of hydrometeors might be also influenced by electric field in the cloud.Low positive Zdr values together with low Zh values correspond to the presence of dry snow, while high Zdr values indicate the presence of pristine ice crystals.
Figure 6f shows that the part of the cloud where W is larger than its surroundings looks like a fountain throwing flakes upward in the central area (thus oriented more vertically) and making the flakes fall and then settle to lower levels further out (thus oriented more horizontally).This can indicate a region of higher turbulence.
The values of Kdp are not available in the region above 2,250 m.The values of Rhohv above 2,500 m are quite low (between 0.8 and 0.9) indicating non-uniform meteorological targets.Due to the low Zh values, predominantly snowflakes or irregular ice crystals can be expected in this region which is in good agreement with that described above.
The measured radar data at lower heights above the radar (up to 2 km) suggest that there is a mixture of different types of hydrometeors including graupel and maybe also hail.Therefore, the resulting environment is suitable for cloud electrification.
These inferences are partly consistent with the official observer report which indicated snow grains falling to the ground.Figure 7 shows the result of the classification of hydrometeors by the XClass algorithm based on the X-band radar data.It indicates that the hydrometeors which occurred mainly near the Milešovka observatory were snow, ice and graupel.This generally agrees with the findings of Figure 6 and the official observer report.
Data measured by the radar MIRA 35c are summarized in Figures 8 and 9. Figure 8 shows the classification of hydrometeors by HClass algorithm for the time interval from 2300 to 2330 UTC.Based on the analysis of the Doppler spectra, the HClass algorithm can identify multiple hydrometeors in a single radar bin, which makes it fundamentally different from the algorithm applied to the FURUNO WR2120 radar data.Therefore, when displaying the classification results, we show the hydrometeors in groups.Specifically, we use the following groups of hydrometeors: CR for cloud or rain droplets, IS for ice or snow, GH for graupel or hail, and ALL for all hydrometeor classes.It should be noted that the classification is done in vertical columns corresponding to different times with a time step of 2 s.The classification (Figure 8) shows that prior to the recorded electrical discharge, the cloud consisted of GH and likely IS and supercooled water.This hydrometeor mixture was identified at heights up to 2,000 m above the radar.Such a combination of hydrometeors is suitable for charge separation by collisions (e.g.Takahashi, 1978;Helsdon et al., 2001;Mansell et al., 2005;Rakov, 2016;Phillips et al., 2020) and the generation of electric field that can be strong enough to cause a discharge.So both the radars suggested presence of a mixture of hydrometeors in the cloud before and during the thunderstorm, which is supposed to be important for cloud electrification, especially in the initial phase.
The results obtained from the analysis of the data measured by the disdrometer (Figure 3) do not fundamentally contradict the classification of hydrometeors derived from the measurements of the MIRA 35c and FURUNO radars.However, although classifications of hydrometeors by the disdrometer and the two radars are consistent, the observer indicated snow grains falling to the ground.
Figure 9 shows the temporal evolution of the vertical structure of the storm for Ze, LDR, Va and Rhohv values from 2300 to 2330 UTC.It is evident from the Ze values that the discharge occurred at the very end of the passage of the storm, that is, after the peak of high reflectivity values.This is consistent with the FURUNO radar data.The LDR confirms the change of near-to-ground temperature caused by the passage of the cold front interface.At 2316 UTC, a melting layer is visible at about 200 m above the ground, indicating a positive temperature at the ground, while the melting layer disappeared when the cold front brought cooler air with temperatures below 0 • C at the ground.This is consistent with Figure 3.However, the temperature drop resulting from the LDR precedes the temperature change measured at the ground by several minutes.This may indicate that the near-surface air mass cooled later than the air above.
The higher LDR values (from −20 to −15 dB) at a height of 1,000 m above the radar indicate that from the bottom point of view, the hydrometeors are more spherical in shape than in the surroundings.This may indicate the existence of graupel or alignment of ice crystals in the electrified field (Caylor and Chandrasekar, 1996;Ryzhkov and Zrni ć, 2007;Hubbert et al., 2014;Melnikov et al., 2018; F I G U R E 7 Classification of hydrometeors using XClass algorithm in the vertical cross-section no. 2 (see Figure 5) on 4 February 2022, at 23:19:16 UTC.The radar position is at (0,0).Classes are the following: DS, dry snow; GR, graupel; HA, hail; IC, ice; LR, light rain; NN, no data; No, hydrometeor cannot be recognized; RR, rain; WS, wet snow.[Colour figure can be viewed at wileyonlinelibrary.com]Sokol et al., 2020).The alignment of ice crystals probably happened at higher elevations where the classification algorithm did not recognize graupel but ice or snow.
As far as we know, behaviour of Rhohv is not well-studied for radars that transmits at only one polarization and receives signal at two polarizations (the case of MIRA 35c).In Figure 9, the Rhohv values close to 1 are at places where LDR is increased.Therefore, we suppose that these values may suggest alignment of ice crystals by the electrified field (Melnikov et al., 2018) at the top of thunderstorm, similarly as in the case of LDR.

General description of the cloud electrification model CEMW
We used our cloud electrification model CEMW (Popová et al., 2022) to simulate electrification of the thundercloud during the winter thunderstorm on 4 February 2022.Our intention was primarily to use the estimated vertical profiles of temperature, humidity, pressure, and horizontal wind components to simulate the evolution of the electric field and to determine whether CEMW gives conditions suitable for a discharge.Apart from the vertical profiles, the CEMW needs an initial impulse (Section 6.2) to form a storm.
Model CEMW makes use of dynamics and two-moment cloud microphysics from the Wisconsin Dynamic and Microphysical Model 2 (WISCDYMM-2) described by Wang et al. (2016).The microphysics considers the following hydrometeor types: cloud droplets, raindrops, cloud ice, snow and graupel.The model assumes a flat terrain and uses a simplified description of the processes taking place very close to the ground.
The CEMW uses open radiative lateral boundaries, Rayleigh relaxation for the upper boundary, and non-slip condition for the bottom boundary for meteorological variables.
CEMW considers seven variables for storm electrification.It takes into account the number of positive and negative ions and the electric charge of the five hydrometeors.CEMW considers the following electrical processes: changes of electric charges of the hydrometeors (included in cloud microphysics), ion equation, induction, and discharges including corona discharges.The model uses boundary conditions, including those for the electric potential, in the same way as Mansell et al. (2005).
CEMW applies Takahashi's scheme (Takahashi, 1978) for change of the charge of hydrometeors due to collisions.The lightning flash scheme used in CEMW is based on that proposed by Barthe et al. (2012) and is similar to the version used by Popová et al. (2022).Details about the model can be found in Popová et al. (2022).

Application of CEMW to the winter thunderstorm
We applied CEMW twice.The two versions differ in the condition on the maximum value of the distributed charge after collision.In the first version denoted integration CEMW, we removed that condition completely, while in the second version called CEMW-lim we used the condition applied by Popová et al. (2022) based on previous work (e.g.Mansell et al., 2005), where limitations for charging were applied to prevent large charging.The maximum magnitude was limited to 50 fC for rebounding graupel-snow collisions per second and 20 fC for graupel-crystal interactions per second.
The lightning flash scheme used in CEMW is based on that proposed by Barthe et al. (2012) and similar to the version used by Popová et al. (2022).Since we wanted to find out if and under what conditions the model will generate a flash, we focused on whether the conditions of our algorithm for a flash are met.The basic condition for the grid point to be able to trigger a flash is that the magnitude where z is the altitude [km] and K = 0.8, a subjectively defined coefficient reflecting uncertainty in calculating E. K includes the impact of limited model resolution on calculated E. The corresponding grid must be within cloud which is defined by total mixing ratio greater than 1.e −7 [kg⋅kg −1 ].
The algorithm triggering lightning differs in two points from the one used in Popová et al. (2022).First, the triggering point has to be in the 3rd and upper vertical levels from the model bottom.We introduced this condition because during the studied thunderstorm, the modelled cloud base reached the ground, which is consistent with observations, but imperfect model boundary conditions at the bottom concerning mainly ions caused artificially higher values of E near the ground.The second more significant change to Popová et al. (2022) is that we did not allow the discharge to extend only vertically in the z direction, as a bidirectional leader does; instead we allowed that the discharge spans also horizontally in x and y directions.This change is motivated by the fact that the charge field of thunderstorms may have a significant horizontal structure (e.g.Brothers et al., 2018;Popová et al., 2022), which agrees with some observations.We ran the model in three domains with different grid resolutions: We integrated the model for up to 35 min.The model domains were set so that the x-axis corresponded to the direction of the front motion (Figure 5).

Initial data and starting up of the model
Initial data, that is, vertical profiles of temperature, humidity, pressure, and horizontal wind velocity components, were used from the forecast of the ICON-EU model (Zängl et al., 2015) valid for 2300 UTC.The model started F I G U R E 9 Measured Ze, linear depolarization ratio (LDR), vertical air velocity (Va), and co-polar correlation coefficient (RHO) by the integrating at 2100 UTC and the vertical profiles were taken from a model grid point located nearest to the Milešovka observatory with a similar elevation.Figure 10 shows the skew T diagram of the initial data.It shows slightly unstable vertical profile of temperature up to 3.5 km approximately where a pronounced temperature inversion is evident, which constrains further vertical extent of clouds.This is consistent with radar and satellite data presented in Section 5.The atmosphere is saturated below the inversion layer.Wind speed is relatively high with prevailing direction from west-northwest.
In the case of modelling summer storms, a warm air bubble is used and placed at the height of the expected cloud base to initiate storm development.Figure 11 shows the temperature initial distribution at the lowest model level.In the present case, the impetus for the convective storm was a cold front, which, along with slightly unstable layer, led to the formation of the storm and induced discharges.When simulating this mechanism (i.e.insertion of cold air beneath the warmer air) by our model, vertical motion was too weak to generate a storm and initiate lightning.Thus, we used, as in the case of summer convective storms, a warm air bubble to initialize the storm but we placed the bubble centre 100 m above the surface.The use of bubble initialization is also supported by the fact that the air was unstable from the ground up to the inversion layer.The bubble radii were 10,000 m in both x and centre varied slightly for each model domain (M1-M3), but this had negligible impact on the results.
Although the warm air bubble initiation does not simulate the cold front conditions at the surface, the character of the model precipitation is qualitatively similar to the measured data both in terms of temporal and spatial evolution.The precipitation field is belt-shaped with a width of about 10 km, which roughly corresponds to the belt width with reflectivity above 20 dBZ (Figure 2).Based on Figure 12, which shows the distribution of precipitation in the model domain, the maximum model precipitation is 1.75 mm⋅hr −1 , which is approximately 60% of the value measured by the rain-gauge at the Milešovka observatory (Figure 3).The main component of precipitation is graupel (about 90%) and partly rain/cloud drops, while snow and ice are not present.The maximum model precipitation duration is 22 min and the observed measured precipitation duration was between 20 and 40 min based on 10 min gauge measurements (Figure 3).
When creating the initial vertical profile, we neglected the wind speed in the y-axis direction (by setting v = 0) in most of the cases, which simplified displaying of the model variables and the model outputs.Although vertical wind shear is known to affect storms, we did not observe any major physical influence on the results.The warm air bubble had a homogeneous structure in the y-axis, therefore the use of v = 0 was logical.Model fields calculated with v ≠ 0 (not shown) were affected by the airflow trajectory making the storm core shift towards the left edge of the model domain in the direction of the airflow.Therefore, most important processes occur near the boundary of the model domain, which could negatively influence the results, and at the same time, v ≠ 0 complicates the presentation of the model results.
Note that we subtracted 15 m⋅s −1 in the x-direction wind component to reduce the size of the model domain.We integrated the model for up to 35 min, but we focused on the evolution of the model fields up to the time when the model conditions triggered a discharge.

Model results for the configuration M1
For presenting the model results, we show the vertical cross-sections passing through the centre of the model domains in the x-axis direction (Figure 11), because the significant processes took place in these areas.
Figures 13 and 14 show resulting charge distribution, charge distribution of hydrometeors, and charge distribution of ions.They also present hydrometeor distribution in the cloud and vertical air velocity in the cloud.Both figures give results based on CEMW run. Figure 13 displays the state after 16 min 40 s of integration, while Figure 14 depicts the state after 33 min 20 s of integration, that is, just before the discharge was triggered in the model.
Both figures evince that the temperature inversion prevents further vertical development of the cloud.Inside the model cloud there are all hydrometeor types, including graupel, ice and snow, which are important for charge separation.Vertical velocity values can be considered reasonable for winter storms.They increased with integration time and reached up to 3 m⋅s −1 at the moment when the of the lightning discharge and in small areas; the rest of the radar velocities were close to 0 m⋅s −1 .
The main difference between the charge fields between Figures 13 and 14 is that the effect of charge separation is obvious in Figure 13.It distributes the charge between graupel on one hand and snow on the other hand.This separation is not visible inside the cloud in Figure 14, where the structure of the charge distribution is very similar for all hydrometeor types.This does not apply to the leeward side of the storm and the ions.
The charge field structure is vertical in the core of the storm where we see the highest vertical velocities.This charge structure is quite different from the smooth idealized fields mainly horizontally oriented, known as dipole or tripole.In our case, a more complicated structure of the charge field is created.At the beginning of the model integration, a charge field has a more or less horizontal structure with three layers differing in sign which is a consequence of the charge separation by collisions.As the storm develops, the vertical wind speed increases in the centre of the storm, which differs significantly from the vertical speed at the sides of the storm and outside of the storm.The resulting motion field affects the shape of the distribution of other model quantities, including those carrying electric charge.Therefore, layers with a predominantly vertical structure corresponding to the vertical velocity structure are formed with different electric charges as a consequence of modelling interactions among hydrometeors and related electrical processes.In contrast, in regions outside of the storm, the model quantities are homogeneous, including the charge fields.
However, if we rotate the field in the direction of the storm axis (vertical velocity field), then the charge fields consist of layers composed of opposite charges and significantly inclined downwards to the right.The structure of the arrays explains why it is appropriate to consider not only the vertical axis in the flash propagation but also the horizontal one.Note that CEMW simulated the first flash after 33 min 28 s of integration at a height of 2,200 m above the surface and the flash leader lay in the horizontal direction.In the case when the limit on the charge rebounding is applied, CEMW-lim did not develop suitable condition for lightning within 35 min of integration.Charge fields of CEMW-lim have similar structure as CEMW but both local extremes and gradients are lower (Figure 15).

Impact of the horizontal resolution on model results
It is well known that the model resolution affects all model quantities.The increase in the horizontal and vertical resolutions is related to the change in the integration time step, which is required to sustain the numerical stability of the model, and should lead to more accurate results.Generally, as the resolution increases, the gradients and local minima and maxima of the modelled quantities increase and the field structures get more pronounced.
In addition to the effect on numerical schemes in our case, the model resolution has a specific impact on the occurrence and evolution of lightning.This is because some processes depend on exceeding the threshold values, which are empirically/subjectively selected and should depend on the model resolution.The Etrig function was derived from the measurements.However, the model operates with discrete values, which typically have lower maxima than continuous arrays of values, and it is reasonable to expect that the maxima of models with higher horizontal and vertical resolution have higher maxima than models with lower resolution.Therefore, the value of K should depend on the model resolution, which we did not apply here, as further research would be needed to select an appropriate K value.
Figure 16 compares total charge distribution and vertical wind velocities for M1, M2 and M3.The displayed fields correspond to the time close to the first modelled discharge.The reason why it is not exactly at the time of the discharge is that the model outputs are always recorded after several model time steps.Figure 16 also shows the vertical velocity which varies from M1 to M3.The basic structure of the wind fields is the same for all three model resolutions, but there are clear differences in maximum values which reach about 3, 3.5 and 4.5 m⋅s −1 for M1, M2 and M3, respectively.The obvious difference is in the detailed structure of the storm with high positive speeds.This area stretches almost vertically for M1, which has the lowest model resolution.For M2 and M3, it is less vertical, and in addition, M3 reaches maximum speeds clearly lower in altitude than M1 and M2.Thus, it is evident from Figure 16 that the model resolution affects the dynamic variables of the model.
The distribution of charge signs in the vertical section is similar for all models, but the values of the charges are clearly different.The M1 needs more time to create conditions suitable for discharge.Here the relationship for calculating Etrig and the fixed K play a role, because it is undoubtedly more difficult to achieve the necessary E values.Differences in vertical velocity fields are also reflected in charge fields.The similarity of charge fields between M2 and M3 is clearly higher than the similarity between M1 and M2, or M3.

Impact of charge separation scheme on model results
The applied algorithm of charge separation by collisions also fundamentally influences the model outputs.For instance, the major influence comes from the limit of allowed amount of charge separation in one model time step or per second, often used to prevent large charging (Figure 15).Figures 17 and 18 show possible values of the amount of charge separation produced by the applied Takahashi (1978) scheme and the impact of the application of charge separation limits (i.e.CEMW vs. CEMW-lim).Figure 17 shows the maximum relative Etrig, which is the ratio of the maximum value of electric field calculated in grid points in the inner model domain using Equation ( 1) in dependence on model lead time for M1, M2 and M3.The maximum relative Etrig is displayed for both CEMW and CEMW-lim runs.The maximum relative Etrig must exceed 1 in a grid point where electric discharge can appear.The inner model domain does not include grid points which are closer than 10 grid points from the model lateral boundaries or are below the third grid point in z direction from the model bottom.The aim of this restriction is to exclude grid points near boundaries where the electric field can be affected by imperfect model boundary conditions.
The results clearly show that limits applied to charge separation slow down the process of cloud electrification.In the case of M1, the applied limit means that no suitable conditions for the formation of the discharge occur within the 35 min of integration, while in the case of M2 and M3 the discharge occurs but with a delay of 7 and 5 min, respectively, as compared to the integration without limits.
The general tendency of the displayed curves in Figure 17 is that they increase with time, but their curves are not monotonous.We believe that the reason is the discretization of the model and the nonlinear processes that affect the Etrig values.It is worth noting that all curves gradually slowly rise up to 25 min of model integration and then sharply increase regardless of the model used.This can be explained by the fact that after 25 min, the convection is strong enough to form a storm.
Figure 17 also shows that there is no close dependence between the maximum relative Etrig and the model resolution (results for M2 and M3 are very similar) and that the dependence can be influenced by the applied limits on charge separation.In our opinion, this behaviour should be attributed to the model numerical algorithms including discretization of the model domain and the nonlinear processes that affect the Etrig values.The results obtained may also be influenced by the particular case-study (we did not record any other winter thunderstorm at the observatory), because this thunderstorm was definitely not typical even in wintertime.Therefore, we are planning to perform similar analyses for our dataset of common summer thunderstorms.
Figure 18 shows maximum magnitudes of charge separation (qsep added to the charge of graupel) per one collision obtained for M1 during the integration time of CEMW.The charge separations graupel-ice and graupel-snow are evaluated separately.Collisions graupel-snow give much greater values than collisions graupel-ice.For both types of collision, the maximum and minimum values as well as mean values of positive and negative charge separation are shown.It is clear that for both types of collision, the maximum and minimum values of the transferred charge are significantly higher (up to two orders of magnitude) than the applied limits.Even the averages are clearly higher than the applied limits (Section 6.1).
The calculation of the qsep according to Takahashi's scheme (1978) depends on temperature and cloud water content and on several parameters/characteristics derived either from the model values (terminal velocities of hydrometeors) or subjectively set for coefficient  (Mansell et al., 2005).Thus qsep may depend on the cloud microphysics used in the model.We plan to look at how the individual components of qsep calculation affect the final value of qsep.To do this, however, we prefer to use data from typical summer thunderstorms rather than the winter thunderstorm studied in this article.

CONCLUSIONS
The present study deals with the analysis and modelling of a winter thunderstorm which was unusual in the Central European context and caused an electrical discharge over the Milešovka meteorological observatory.Lightning network EUCLID registered two strokes immediately following one another at 2320 UTC on 4 February 2022.Those two strokes were recorded at a distance of less than 200 m from the observatory.One of them was also reported by an observer having position at the station.
The study consists of two parts.In the first part, data from the Meteosat Second Generation satellite, standard meteorological measurements at the station, forecasts by the numerical prediction model ICON-EU, and data from the X-band Doppler polarimetric radar and the Ka-band Doppler polarimetric vertical profiler were used to analyse the state of the atmosphere at a time of lightning.The second part of the study focuses on simulating the evolution of the winter thunderstorm with an emphasis on cloud electrification and suitable conditions for lightning.
The thunderstorm occurred on a cold front crossing the Milešovka from the northwest.This cold front passage is clearly visible in both the station and the radar data.Measurements of the Meteosat Second Generation satellite and of both radars show that the cloud height was quite low for a thunderstorm (about 3.5 km above sea level).According to the data of the ICON-EU numerical weather prediction model, there was a temperature inversion at this altitude preventing further vertical development of cloud.
Radar data are useful for the identification of the storm's internal structure and its evolution before and after lightning.Our results based on radar data can be summarized as follows: • The discharge occurred at the very end of the passage of the storm core showing high reflectivity values over the Milešovka observatory.We suppose that the reason is that the process of cloud electrification requires some time to produce a sufficiently strong electric field that leads to lightning.
• Classification of hydrometeors for both radars confirmed that the hydrometeors, which are considered necessary for charge separation by collisions (graupel, cloud or rain water, ice or snow), were present in the thundercloud.This mixture of hydrometeors was found in the data measured by the cloud vertical profiler from the ground level up to 1,000 m for a long time before the discharge occurred.The vertical extent of this mixture increased sharply up to 2,000 m 5 min prior to the discharge, approximately.On the contrary, that mixture disappeared completely after the discharge.
• Radar data showed inhomogeneity in values of the majority of measured/derived quantities at altitudes above approximately 2,000 m about 5 min before the discharge.The cloud profiler showed short-term increase in updraught velocities at different altitudes.At the same time, an increase in LDR and Rhohv was visible at higher altitudes.The causes of it cannot be clearly determined; however, a possible explanation is that the created electric field may affect the organization of ice particles, which was suggested previously in literature.
X-band radar cross-sections indicated large variability of radial velocity of targets and probably rotational air movements at 2,000 m and above the radar in the vicinity of the Milešovka observatory.
Model results are burdened with some uncertainty coming from the way the storm was initiated.We used a method often used for summer thunderstorms, that is, the warm air bubble.The reason for this was that the cold front passage could not be well modelled with our model.Unlike summer storms, we placed the centre of the warm air bubble just above the surface (100 m).In our opinion, while this initialization does not realistically model the cold front passage close to the Earth's surface, the model behaviour at higher altitudes is realistic.
Although the main focus of this study was to show whether the model is capable of producing lightning under given meteorological conditions, the model-simulated storm development has some similarities with the general knowledge about winter thunderstorms presented in the Introduction (Section 1).Lightning appears when the storm is in its mature stage characterized by maximum vertical velocities and mainly graupel falling to the ground (the model does not consider hail).This state occurs for approximately 18 min, which approximately corresponds to 15 min indicated by Altaratz et al. (2021).
Results for model versions M1 (lowest resolution), M2, and M3 (highest resolution) can be summarized as follows: • The basic structure of model fields is similar for M1, M2 and M3 but there are differences in local extremes depending on the resolution.It was found that the model resolution affects the storm dynamics, which was apparent in vertical air velocity fields.
• All three model configurations gave similar structure of charges.Bands with opposite charge had a distinct vertical component.When the axis was taken to be aligned with the storm axis, then the bands were inclined for M2 and M3.This made us modify the lightning algorithm to consider not only vertical propagation of the leader but also horizontal propagations in x and y directions.
• All the model configurations generated lightning within 35 min of integration time when the applied condition on the magnitude of charge separation was removed.
The model with the lowest resolution M1 triggered a discharge after 30 min 50 s, with M2 and M3 about 5 min sooner.When the restriction on the magnitude of charge separation was applied, M3 and M2 produced lightning 5 and 7 min later, whereas M1 did not develop any lightning within the 35 min of integration.
• Takahashi's scheme gave significantly larger values of charge transferred during collisions of graupel with ice and snow than the commonly used maximum values.We found that, in absolute value, positive changes of charges are greater than negative contributions for both types of collisions.We got that the maximum absolute values of transferred charges in single model time steps may be about two orders of magnitude larger for graupel-snow collision and much smaller for graupel-ice collision.However, these maxima are reached at individual grid points.Average values are much smaller but still several times larger than the usual limits.The introduction of limits clearly slows down the process of discharge generation, more significantly in the case of lower model resolution.This is logical because lower resolution leads to smoother model fields.We believe that if such a constraint on the transmitted charge has to be introduced, then it should depend on the temporal and the spatial resolution of the model.
The calculation of transferred charge based on Takahashi's scheme depends on several components including the cloud microphysics used in the model.In the future, we plan to look at their impact on the resulting electric field by using a dataset of observed summer thunderstorms, and will present the results in the upcoming study.

F
Meteorological observatory Milešovka.[Colour figure can be viewed at wileyonlinelibrary.com] that winter thunderstorm.Conclusions and future plans are given in Section 7.

F
I G U R E 2 Radar fields at constant altitude plan position indicator (CAPPI) 2 km ASL during the winter thunderstorm: (a,b) Zh, (c) RhoHV and (d) Zdr.The time of the beginning of plan position indicator (PPI) scans is shown in the title in the format YYMMDD for year, month and day, respectively, and HHMMSS for hour, minute and second, respectively, UTC.The black point in the centre of the panels indicates the position of the Milešovka observatory and the axes show the distance from the observatory.The black line represents the boundary between the Czech Republic and Germany.Horizontal and vertical axes show distances from Milešovka in m.The dashed arrow shows the approximate movement of the cold front.[Colour figure can be viewed at wileyonlinelibrary.com]

F
Values of brightness temperature [K] in the infrared channel IR 10.8 μm measured by Meteosat second generation (MSG) satellite on 4 February 2022, at 2320 UTC.The black cross depicts the position of the Milešovka observatory and the black line represents the border of the Czech Republic.[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 5 Vertical cross-sections of Zh at 23:19:16 UTC (RHI scans).Red arrows in the upper-right panel show the orientations of the cross-sections, while the green arrow in the upper-right panel shows the approximate movement of the cold front (i.e.cloud belt) at constant altitude plan position indicator (CAPPI) 2 km.The direction of the cold front corresponds to the wind direction measured at the surface at the time of the recorded lightning.The radar position is at (0,0).[Colour figure can be viewed at wileyonlinelibrary.com]

F
Classification of hydrometeors in vertical columns of the MIRA 35c radar using HClass algorithm on 4 February 2022, from 2300 to 2330 UTC.Classes are the following: ALL, all hydrometeors; CR, cloud or rain droplets; GH, graupel or hail; IS, ice or snow.The time step is 2 s.The vertical dashed line denotes the time of the recorded lightning.[Colour figure can be viewed at wileyonlinelibrary.com] of grid electrical field E exceeds a threshold value Etrig [kV⋅m −1 ], Model configuration M1: 128 × 78 × 107 grid points in x, y and z directions, respectively, with mesh sizes of dx = dy = 500 m, and dz = 200 m, and a time step of dt = 2 s. • (ii) Model configuration M2: 128 × 78 × 210 grid points in x, y and z directions, respectively, with mesh sizes of dx = dy = 250 m, and dz = 100 m, and a time step of dt = 1 s.• (iii) Model configuration M3: 128 × 78 × 305 grid points in x, y and z directions, respectively, with mesh sizes of dx = dy = 150 m, and dz = 50 m, and a time step of dt = 0.5 s.
MIRA 35c radar from 2300 to 2330 UTC.The vertical dashed line shows the time of the recorded discharge.[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 10 Skew T diagram showing model initial data derived from the ICON-EU data valid at 2300 UTC. Green and red curves show vertical profiles of temperature and dew-point temperature, respectively.[Colour figure can be viewed at wileyonlinelibrary.com] z direction and 90,000 m in the y direction.The centre of the bubble was 2 K warmer than the original temperature.The shape of the bubble emulated the band shape of the observed thunderstorm.The x-axis location of the bubble F I G U R E 11 Temperature initial distribution at the lowest model level.The x-axis corresponds to the flow direction and the dashed line shows the position of the vertical cross-section used for visualization of the model outputs.

F
I G U R E 12 Distribution of total precipitation [mm⋅hr −1 ] calculated by CEMW (M1 configuration) at the central part of the model domain.The x-axis corresponds to the flow direction and the dashed line shows the position of the vertical cross-section used for visualization of the model outputs.

F
I G U R E 13 Vertical cross-section of CEMW: (a) total charge distribution, (b) hydrometeor distribution in the cloud and isotherms 0 and −15 • C, (c) vertical velocity in the cloud oriented upwards, (d-j) charge distribution of hydrometeors, and (k) charge distribution of positive and negative ions.The forecast time is shown in the MM:SS format.In panels (a) and (d-k), solid and dashed black contours show positive and negative vertical velocity, respectively.Contour values are −2, −1.5, −1, −0.5, −0.1, 0.1, 0.5, 1, 1.5 and 2 m⋅s −1 similar, as in panel (c).In panel (b), a hydrometeor is shown at a grid point if its mixing ratio exceeds 1.e −6 [kg⋅kg −1 ]. [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 14 The same as in Figure 13 but for the integration time of 33 min 20 s. [Colour figure can be viewed at wileyonlinelibrary.com] lightning discharge was triggered.Compared to the estimated maximum of vertical velocity derived from the radar data (Figure 9), the model maxima are much lower (less than a half), though the comparison is difficult because the spatial resolution is different and the radar measured higher speeds in very short periods of time around the time F I G U R E 15 Total charge fields in nC after 33 min 20 s of integration by (a) CEMW and (b) CEMW-lim model runs.Solid and dashed black contours show positive and negative vertical velocity oriented upwards, respectively.Contour values are −2, −1.5, −1, −0.5, −0.1, 0.1, 0.5, 1, 1.5 and 2 m⋅s −1 .[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 16 Total charge q for (a) M1, (b) M2 and (c) M3, and vertical velocities for (e) M1, (f) M2 and (g) M3.The displayed fields correspond to the time close to the first discharge developed by the model runs.Solid and dashed black contours in panels (a-c) show positive and negative vertical velocity oriented upwards, respectively.Contour values are −2, −1.5, −1, −0.5, −0.1, 0.1, 0.5, 1, 1.5 and 2 m⋅s −1 .[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 17 Maximum relative Etrig (q) as the ratio of the maximum value of the electric field calculated in grid points in the inner model domain using Equation (1) as dependent on model lead time for M1, M2 and M3.The solid and dashed lines show results of CEMW and CEMW-lim model runs, respectively.[Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 18 Charge separation values (qsep added to the charge of graupel) per one collision obtained for M1 as dependent on the integration time for graupel-ice and graupel-snow collisions.Solid coloured lines show the maximum and the minimum values for the two types of collisions, while dashed lines show mean values of positive and negative charge separations.The two horizontal black lines indicate 50 and −50 fC values.[Colour figure can be viewed at wileyonlinelibrary.com] Jana Popová conceptualization, writing, review & editing, validation.Zbynek Sokol conceptualization, methodology, investigation, writing, software, visualization.Pao Wang methodology, writing, review & editing.Jaroslav Svoboda investigation, software.