1 INTRODUCTION

Winter thunderstorms are rare in Central Europe and some other parts of the world such as China. They are so rare and were regarded as a bad omen that it was required to be written in records in official chronicles of Ancient China whenever they occurred (Wang, 1980). Although they are usually not related to significant damages as compared to summer thunderstorms, from a meteorological point of view they are an interesting phenomenon worth investigating. Most attention in storm research is, as expected, paid to summer storms, though 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 5 deals with lightning data, while Section 6 presents basic meteorological data. Sections 7 and 8 show our results. Section 7 describes the vertical structure of the winter thunderstorm derived from data of the two radars, whereas Section 8 focuses on the application of our cloud electrification model to simulate that winter thunderstorm. Conclusions and future plans are given in Section 15.

2 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.

3 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.

4 METEOROLOGICAL CONDITIONS OF THE WINTER THUNDERSTORM ON 4 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.

5 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; 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.

6 MODELLING WINTER THUNDERSTORMS USING A CLOUD ELECTRIFICATION MODEL

6.3 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 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 7. 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 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 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⁻¹, 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⁻¹ 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.

6.4 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⁻¹ at the moment when the 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 of the lightning discharge and in small areas; the rest of the radar velocities were close to 0 m·s⁻¹.

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).

6.5 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⁻¹ 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.

6.6 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 8).

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.

7 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.

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).

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.

AUTHOR CONTRIBUTIONS

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.