Radar-based assessment of hail frequency in Europe

In this study we present a unique 10-year climatology of severe convective storm tracks for a larger European area covering Germany, France, Belgium, and Luxembourg. For the period 2005-2014, a high-resolution hail potential composite of 1×1 km is produced from two-dimensional radar reflectivity and lightning data. Individual hailstorm tracks as well as their physical properties, such as radar reflectivity along the tracks, were reconstructed for the entire time period using the Convective Cell 5 Tracking Algorithm (CCTA2D). A sea-to-continent gradient in the number of hail days per year is found to be present over the whole domain. In addition, the highest number of severe storms is found on the leeward side of low mountain ranges such as the Massif Central in France and the Swabian Jura in Southwest Germany. A latitude shift in the hail peak month is observed between the northern part of Germany where hail occurs most frequently in August, and southern France where the maximum of hail is two months earlier. 10 The longest footprints with high reflectivity values occurred on 9 June 2014 and on 28 July 2013 with lengths reaching up to 500 kilometers. Both events were associated with hailstones measuring up to 10 cm which caused damage in excess of e2 billion.

criterion on single-polarisation radar data to retrieve hail signals for the period 2002 until 2011. The authors found that hail 60 occurred mostly during May, June and July during the afternoon throughout the country. Despite the use of improved radarbased techniques, most of the studies cited above were restricted to smaller regions or a single country. While some authors have estimated hail frequency from other sources such as Overshooting Tops (an indicator of strong convective updrafts) in satellite imagery (Punge et al., 2014), model data (Mohr et al., 2015b, a;Rädler et al., 2018) or a combination thereof (Punge et al., 2017), radar proxies have a key benefit over overshooting tops as they represent a more direct measure of hail within a 65 storm (presence of large reflectivities). Numerical models, such as weather forecast models or regional climate models (RCM), on the other hand, are not able to reliably reproduce hail due to a high degree of uncertainty in the initial conditions, a lack of knowledge in cloud microphysics, and the high computer costs when running a two-or three-moments microphysics scheme.
The big advantage of satellite or model-based hail proxies is that they cover large areas more or less homogeneously, whereas radar-based climatologies are typically limited to a single country or region. 70 Our study is the first to combine radar observations from multiple countries. Indeed, the objective of our study is to analyze the spatiotemporal variability of hail signals over a 10-year period (2005 to 2014) covering the four European countries of France, Germany, Belgium and Luxembourg. Hail signals were estimated from 2D radar reflectivity available for each country, which permits a homogeneous hail analysis. The results help to identify regions frequently affected by hail and allow us to relate hail frequency to topographic features such as terrain height or the proximity to the sea. A thorough study of hail 75 events gives further insights into the relation between orography and deep moist convection. Improved understanding of these mechanisms and processes is crucial to improve the nowcasting and forecasting skill of hail storms. Finally, as hail constitutes a considerable risk for the insurance industry, improved knowledge about hail frequency and hailstorm characteristics will help to better understand the related risks.
The paper is structured as follows: Section 2 gives an overview of the remote-sensing and reanalysis datasets used for this 80 study. Section 3, describes the combination of radar data with lightning data and the application of the tracking algorithm.
The remote-sensing output is then used to generate European composites at 5-minute time steps. Section 4 assesses the hail variability between 2005 to 2014 in relation to the distance to the sea and the presence of orography near hailstorms. This section also presents results on seasonal and diurnal variations in hail frequency and provides some characteristics of the hail cells. Concluding remarks follow in Section 5.

French radar data
The French radar network operated by Météo-France and the derived radar products evolved constantly through times, mainly 100 via national projects. A brief overview of the French radar network is given hereafter. In 2001, 19 radars were constituting the French radar network. One year later, in 2002, five new radars were added to the network and some of the 19 radars were replaced by dual polarization radars (Tabary et al., 2006, Bousquet et al., 2008. In 2005, 24 radars were in operation including 19 C-band radars and five S-band radars (both with a radius of up to 120 km). Two years later, in 2007, the radar stations of Toulouse in southwestern France and Trappes (near Paris) were renewed (Tabary, 2007) but this replacement did not affect 105 the radar national composite. During the period from 2007 to 2011, the radars of Plabennec located in northwestern France, Abbeville in northern France, Nimes in southern France and Grèzes in the southern part of central France were replaced as well with dual polarization radars. In 2014, five X-band radars with an average coverage radius of 50 km were added to the French national radar composite in the Alpine region (Beck andBousquet, 2013, Champeaux et al., 2011). As the data from the X-band radars were only recently implemented into the French national composite (Yu et al., 2018), only S-and C-band radars 110 were considered in this study. Note that the radar stations of Avesnois (located in northern France) and Réhicourt-la-petite in Lorraine (labelled as 7 in Figure 1) cover a large part of eastern France, and permit to integrate Luxembourg completely, as well as a significant part of Belgium, into the French national composite. Concerning the scanning strategy, four to six scans are performed every 15 minutes at elevation angles ranging from 0.4°u p to 15° (Figueras i Ventura and Tabary, 2013). Only the lower elevation angles below 2.7°are scanned every 5 minutes.

115
The spatial resolution of the composite is 1 × 1 km 2 with a size of 1536 × 1536 grid points for each image referred to a plane Cartesian coordinate system (Tabary et al., 2006). Radar data from all stations are pre-processed via an algorithm (Tabary, 2007) named Castor2 (Figueras i Ventura et al., 2012), which corrects for several errors, such as antenna positioning errors, and quantifies horizontal reflectivity Z h in polar coordinates (Tabary, 2007). During the pre-processing stage, each radar pixel receives a weighted quality index (QI) ranging from 0 to 1 (Tabary, 2007), updated throughout the whole pre-processing chain.

120
The first pre-processing step is to eliminate ground clutter, i.e., fixed echoes at the surface, using Doppler velocity . Then an orographic mask is applied at each elevation angle in order to assess the beam occultation rate. After that, an "anthropogenic" mask, including buildings, trees, or other fixed objects in the vicinity of the radar is computed with the help of long-term accumulated radar products. These masks allow to remove radar pixels with artificially high reflectivity at each elevation angle. The underestimation of reflectivity above the 0 • C isotherm is taken into account  using 125 vertical reflectivity profiles. Attenuation by oxygen is corrected depending on the wavelength, the elevation, and the distance to the radar (Doviak and Zrnić, 2006). For example, a correction of 1.79 dBZ at 100 km away from the radar site is applied to C-band radars for an elevation angle of 0.4° . One of the last pre-processing step is the correction of the bright band with the help of vertical reflectivity profiles. After performing all the steps described above, individual plan position indicators (PPIs) are combined to 2D composites produced every 5 minutes available for each radar station. All individual radar 130 products are then combined into a national mosaic . For areas with overlapping radar coverage, weighted reflectivity data are computed depending on the distance to the nearest radar . At the boarders of France, radar data from other national weather services are integrated into the French national mosaic. Reflectivity values from the French radar mosaic used in this study were coded in a table and stored in GeoTIFF format, i.e., georeferenced TIFF images.
The resolution is 2 × 2 km 2 from 2004 until mid-June 2009, and a finer resolution of 1 × 1 km 2 is available from mid-June 135 2009 to 2014. For data homogenization, each of the 2 × 2 km 2 composite was interpolated linearly from 2004 to June 2009 to the finer grid of 1 × 1 km 2 .

German radar data
The German radar composites for the period from 2005 to 2014 are provided by the German Weather Service (DWD), which operated a network of 17 C-band radar systems in 2014. During the investigated period, a new radar in Memmingen, southern 140 Germany, was added to the network in 2012 (Puskeiler, 2013). As the horizontal range detection for each radar is 180 km and a maximum distance of 200 km separates the radar stations, an extensive overlap of the detection areas permits almost a complete coverage of the German territory. Only some peripheral regions are not well covered by the composite, for example, in the far North near the Danish border or in southeastern Bavaria ( Figure 2). In the complex terrain of southern Germany, weather radars are preferably located on hills and mountains to minimize beam shielding by orography. Concerning the scanning strategy, the 145 lowest elevation angles between 0.5°and 1.8° (Bartels et al., 2004) at that time were scanned every 5 minutes, whereas a complete volume scan took 15 minutes. The maximum reflectivity values of the lowest elevations are used for the national 2D reflectivity composite. The steps of the pre-processing are similar to those of Météo-France, and include among others: An elimination of clutter pixels using a clutter filter, orographic shading correction using an elevation model. After the preprocessing, the local radar data are merged into the German national composite. For areas with overlapping radar coverage 150 the maximum reflectivity value from all radar scans is used in the composites, while for the neighboring regions of foreign countries, a weighted adjustment is performed between radar products from other national weather services and the German rain-gauge dataset (Kreklow et al., 2020). The quality of German radar data has improved over the last decade with continuous algorithm corrections and adjustments (Kreklow et al., 2020) used for RADOLAN (Radar-Online-Aneichung, which means Radar-Online Adjustment) and can be assessed by a quality flag provided for each pixel on the reflectivity product. The spatio-155 temporal resolution as well as the time period available for the German radar data is the same as the French, namely 1 × 1 km 2 with a 5-minute time step for the composite and available from 2005 until 2014, so that both data sets can be merged. The encryption of each scan of the DWD network entails near ground reflectivity values (named RX product) in so-called RVP6 units. The advantages of this data are the high temporal and spatial resolution which enables us to properly identify footprints of SCS. RX data are projected on a Cartesian grid so that each grid box is equidistant at 1.0 km. In the end, the German radar 160 composite has a size of 900 × 900 km 2 covering the whole of Germany.

Uniform Pan-European Grid
It is important to note some limitations in both the German and French national composites. Long-term QPE maps for the French national composite reveal some regions with low data accuracy. This is mainly the case for the central and eastern part of the Pyrenees mountains and the entire Alpine region . In the other parts of France, the QI is 165 mostly higher than 90% with especially high QI close to the radar site . Radar data failure, for example during radar calibration, or radar replacement, were estimated by Puskeiler et al. (2016) to be approximately 4.5 ± 3.9% on average (mean ± standard deviation for the German national composite). Furthermore, the combination of the German and French national composites, each calibrated and pre-processed in different ways, may lead to inhomogeinities in relative hail frequency in some regions. Based on manual investigation of several cases with severe hailstorms in the border region 170 between Germany and France, it was found that the signal of the French mosaic is between 0.5 and 1 dBZ lower compared to that obtained from the DWD composite (Schmidberger 2020, personal communication). This uncertainty is acceptable when projecting the two national composites onto a uniform Pan-European Grid. Radar reflectivity data and thus radar-derived hail signals were projected on the same uniform European grid (not shown) with a resolution identical to that of the national radar network (1 × 1 km 2 ). We used the geographic coordinate system WGS84 for the Pan-European Grid and a Lambert Conformal 175 Conic Projection, as recommended by Gregg and Tannehill (1937) and Varga (1990). In the center at about 47 • N and 6 • E, the meridional grid spacing is equal to the zonal direction to minimize the grid distortion.

Lightning data
To remove artificial clutter still present in the data, we additionally implemented a filter based on lightning data, which was already used by Puskeiler et al. (2016). Here we used only cloud-to-ground (CG) lightning (strokes) from the low-frequency 180 lightning detection system BLIDS (BLitz InformationsDienst Siemens), which is part of the EUCLID (EUropean Cooperation for LIghtning Detection) network. The detection efficiency of the system is 96% for strokes with a peak current of at least 2 kA (Schulz et al., 2016). Because the sensors and the algorithm implemented until 2015 had a significantly lower detection efficiency of intra-cloud and cloud-to-cloud lightning according to Pohjola and Mäkelä (2013), these types of lightning were not considered.

ERA5 reanalysis
To assess the mean wind flow during hail days, we used the ERA5 global reanalysis (Hersbach et al., 2020). ERA5 is a new global atmospheric reanalysis recently released by the ECMWF and aims to replace ERA-Interim reanalysis (Dee et al., 2011) whose data extend from 1979 to 2019. For the moment, ERA5 is available from 1979 onwards and will be soon extended 190 to 1950. The ERA5 4D-Var analysis dataset is assimilated by the Integrated Forecasting System (IFS) and is available on a horizontal resolution of 0.25 • on 137 vertical levels every hour.

Correction of erroneous signals 195
Concerning the homogenization of the French and German national composites, several corrections had already been performed by both national meteorological services. Radar reflectivity data still contains noise and systematic errors that have to be eliminated using various approaches. Errors mostly concern individual radar pixels with significantly higher reflectivity values (e.g., more than 70 dBZ) compared to the surroundings. To avoid this problem, reflectivities below 35 dBZ or above 70 dBZ were set as missing values. Following Puskeiler et al. (2016), an additional verification and correction filter was applied for 200 reflectivity values of Z > 45 dBZ with a difference of ∆Z > 5 dBZ to the adjacent pixels. The affected pixel is set to the mean value of its 8 surrounding pixels and this filter was applied to all consecutive radar scans: In addition, if a reflectivity value at a given grid point is at least twice as high compared to the 8 neighboring values and was not present in the scan before or afterward, the reflectivity value is considered an artifact and set to zero.

Lightning filter
Although the radar tracking routine (see next paragraph) includes a clutter filter, several erroneous signals are still present in the radar data. For example, isolated non-meteorological targets such as electronic signals or reflectivities from wind turbines can emerge in radar scans (Steiner and Smith, 2002).
Since hail occurs only in association with thunderstorms (Baughman and Fuquay, 1970;Changnon, 1999;Wapler, 2017), 210 lightning is expected near high reflectivity cores. In addition to the gradient filter described above, we used lightning detections to further remove artificial clutter. If high reflectivity values (Z 55 dBZ) occur during a 24-hour period without lightning, the values at the affected grid points are set to zero. A maximum distance of 10 km was chosen between a lightning discharge location and the pixels with high reflectivity. Distances of 5, 15, and 20 km were also tested; a distance of 5 km led to the disruption of several hail tracks due to gaps in reflectivity values; the other two thresholds affected the results only marginally.

215
An example of the lightning filter application during a hailstorm can be found in Fluck (2018)  The object-based Convective Cell Tracking Algorithm (CCTA2D) permits the reconstruction of tracks of individual convective cells using 2D radar data. The algorithm is based on the tracking algorithm TRACE3D (Handwerker, 2002), originally devel-220 oped and optimized for 3D radar reflectivity from a single radar in spherical coordinates. TRACE3D was further extended to radar reflectivity data in Cartesian coordinates such as those provided by the DWD radar network (Puskeiler et al., 2016). A second version was adapted to 2D terrain-following near ground reflectivity (CAPPI) using both the RX product from DWD and the French mosaic including France, Belgium, and Luxembourg (Fluck, 2018).
The first step of CCTA2D is to identify, regions of intense precipitation (ROIP) delimited by Z ≥ 35 dBZ and to determine 225 the corresponding maximum reflectivity values (Z max ). In order to distinguish individual reflectivity cores (RCs) within each ROIP, a value of ∆Z = 10 dBZ is subtracted from Z max to set the minimum threshold necessary to delimit an single RC (Z rc ).
Thus, the value of Z rc remains the same for all identified RCs inside a ROIP. If Z rc is less than 55 dBZ, the RC is rejected and not tracked by CCTA2D. Two additional conditions are required for an RC to be classified as a potential convective cell and to be tracked by CCTA2D: A minimal area of 5 km 2 is needed to define an RC with at least 3 pixels (km 2 ) of Z ≥ 55 dBZ. The 230 thresholds detailed above to identify potential convective cells in CCTA2D are summarized in Table 1. The 55 dBZ threshold is referred to as the hail criterion according to Mason (1971), and was successfully used in several studies (e.g., Hohl, 2001,Hohl et al., 2002,Kunz and Kugel, 2015. Schuster et al. (2005), for example, found the 55 dBZ to be a good indicator for damaging hail on the ground in Eastern Australia. Puskeiler et al. (2016) estimated a slightly higher threshold of 56 dBZ best separating between days with and without insured damage to buildings, but confirmed the 55 dBZ to estimate at best insured damage to 235 crops. Categorical verification using insurance loss data over a 7-year period in southwest Germany for this threshold yields a Heidke Skill Score HSS of 0.6, a quite high value confirming the detection skill (it should be noted that this value increases to HSS = 0.71 when using an adjusted version of the Waldvogel et al. (1979) criterion requiring 3D radar data). In the same study, Puskeiler et al. (2016) found that the Probability of Detection (POD) reached 0.65 and the False Alarm Ratio (FAR) was 0.4, indicating that 35% of the observed hail events are missing while 40% of those predicted events are false alarms.

240
The second step of CCTA2D is the temporal and spatial tracking of all detected convective cells. The algorithm associates RCs between consecutive radar composites according to the estimated propagation velocity and the position of the RC. Prerequisite of the tracking is that RCs intensity and size, from one time step to the next must exist within a certain search radius for accurate RC assignment and tracking. The search radius is given by the estimated distance of an initial RC displaced during a time step of 5 minutes multiplied by a velocity factor of 0.6.

245
Special attention is given to cell splitting and merging. Cell splitting is a prominent feature of supercells associated with vertical pressure disturbances. In the northern hemisphere, where the hodograph is usually right-curved, right-moving storms tend to be favored compared to left-moving storms. Cell splitting may also occur due to changes in storm intensity that cause a single RC to break up (or vice versa in the case of cell merging). In order to track both cells after they have split, a splitting (merging) option in the tracking algorithm is necessary. Furthermore, without splitting or merging options, the physical 250 characteristics of SCS tracks such as their length or their angle of orientation could be incorrectly computed by CCTA2D. To Minimum number of pixels inside an RC 3 km 2 detect cell splitting, the initial cell (e.g., the "parent" cell) is first spatially displaced to the position of the following cell (e.g., the successor, or "child" cell), and their respective areas are compared (Handwerker, 2002). A split is defined when a cell at time t can be associated with two cells at time t+ ∆t. In this case the largest "child" cell inherits the history of the "parent" cell.
Similarly, a merger is defined when two cells at time t are associated with a single cell at time t+ ∆t. In this case, the largest 255 "parent" cell assigns a history to the "child" cell. The maximum distance between two RC centers that could merge is set to 10 km. Initial and successor areas are then compared, and the successor is placed at the weighted center of all initially detected cores. Merging occurs when the successor area is larger than the initial RC. To avoid reflectivity core crossings or overlapping, each RC is enumerated and recorded separately.
After the construction of entire cell tracks, the composite of maximum reflectivity on a given day does not provide a smooth 260 result, but a rather scattered product. This effect is most pronounced when the cells propagate further than their horizontal extent during a measuring interval. The faster the storms move, the more scattered is the maximum reflectivity projected on a 2D plane. This can substantially reduce reflectivity values between two scans, even though a high-intensity storm crossed the area. A gap of reflectivity values can also appear on radar scans in regions with overlapping radar data, especially on neighboring countries such as in the Rhine Valley. To consider this effect, an advection correction was performed following 265 Puskeiler et al. (2016). A translation of the reflectivity cores is computed from one time step to the next considering the horizontal wind field estimated by CCTA2D along a track. The field of motion vectors are computed and projected on the German and French grid. Each point along a track includes a velocity shift-vector in north-to-south dv and west-to-east du directions. The so-called shift-vector U is denoted as: A method is applied to obtain smoothed reflectivity values along the tracks detected by the algorithm CCTA2D. First, a cluster or "cell" of reflectivity values is detected along the track at a time t. This "cell" includes the maximum reflectivity value at time t as well as its 40 km by 40 km 2 surroundings reflectivity values. After that, the cell is displaced along the track with a radius of 3 km (Puskeiler, 2013). Then, the first "cell" of reflectivity values is shifted forward in time and the second "cell" of 275 reflectivity values is shifted backward in time.  hail criterion of Z ≥ 55 dBZ is fulfilled on a specific day at a single grid point, this grid point is set to 1, otherwise it is counted as zero. The total of all days with hail over the entire 10-year period divided by the number of years yields the radar-based "hail climatology". In accordance with other hail frequency analyses (e.g., Puskeiler, 2013, Nisi et al., 2016, Junghänel et al., 2016, Nisi et al., 2018, the term climatology is used here, even though our investigation refers to a period far below a climatological time scale of ≥30 years. Note that this climatology represents the spatial distribution of convective cells with high reflectivity,

295
but not directly of hail as the 55 dBZ threshold does not guarantee hail on the ground. Similarly, the absence of high reflectivity does not ensure that hail did not occur (See Section 3.3). The term hail days used in the following parts of this study refers to the exceedance of reflectivity, but not to confirmed hail observations.
As can be seen in Figure 3, the spatial variability of hail days per year is very large, but some patterns with distinct minima or maxima can be identified. A close-up investigation of the hail hot spots is detailed hereafter. Figure 4 represents the location of the mean hail days per Sea with a southerly direction, thus impinging the southern and southeastern mountains of the Massif Central at a sharp angle.
Another general westerly flow reaches the western part of the Massif Central. Interestingly, it seems that not only the location of the Massif Central is responsible for the increased number of hail days downstream, but also the flow convergence where the westerly flow meets with the flow coming from the Mediterranean. One may speculate that even without the Massif Central hail days might be increased in that area of low-level flow convergence. The large valleys on the western side of the Massif 320 Central, oriented from southwest to northeast, facilitate the passage of the flow coming from the southwest into the Livradois region. This region with an average number of 3.2 hail days per year is located in an area where the wind vectors converge both in the direction and velocity. In order to better understand the flow characteristics over the Massif Central shown in Figure 4, we calculated the Froude number on radar-derived hail days from ERA5 (Queney, 1948;Smith, 1979) for a region covering the Massif Central entirely and ranging from 44.0°to 46.5°N and from 2.0 to 4.7°E. The Froude number is calculated as follows: where U represents the wind speed perpendicular to the mountain and was computed by applying a density weighted integration over the lowest 2000 m. H is a characteristic mountain height set to 1300 m for the Massif Central region and N is the Brunt-Väisälä frequency. The Brunt-Väisälä frequency N is defined as : Where g is the gravitational acceleration equal to 9.8 m.s −1 , θ v is the virtual potential temperature, and ∂θv ∂z represents the vertical gradient of the virtual potential temperature. In our analysis, we considered the root-mean-square of the Brunt-Väisälä frequency N in order to exclude imaginary values.
The mean Froude number on hail days over the Massif Central from 2005 to 2014 is F r = 0.39 ± 0.3. According to Smith (1979) and Smolarkiewicz and Rotunno (1989), a Froude number below 1 suggest a flow that goes around the mountain rather 335 than directly over it. Thus, it can be assumed that the flow around the Massif Central is deviated by the mountains peaks leading to convergence downstream at low levels on the leeward side of the Massif Central, where the hail hot spot is located. This deflection of the flow is unlikely to show up in 4 due to the fairly coarse resolution of the ERA5 reanalysis.
Several authors have found an increased hail frequency rather downstream than upstream or directly above the mountains. This is for example the case in the Pyrenean region (Vinet, 2001, Berthet et al., 2011, Hermida et al., 2013 2013), near the Black Forest in Germany (Kunz andPuskeiler, 2010, Puskeiler et al., 2016), or in the vicinity of the Alps (Eccel et al., 2012, Nisi et al., 2018. By referring to the studies of Mass (1981)  14 355 Another hail hotspot in the northeastern part of France is found along the northern ridge of the Jura Mountains (labelled E in Figure 1) in Franche-Comté with 2.5 hail days per year.
Note that the Jura mountains represent a natural obstacle frequently triggering thunderstorms (Piper and Kunz, 2017) and hailstorms by orographical lifting (e.g. Langhans et al., 2013, Schemm et al., 2016, Nisi et al., 2018. The Rhone-Alpes (region 11 in Figure 1) is a region likewise frequently affected by hail. This region contains the large 360 Rhone valley and is bordered by the Massif Central in the west and by the Alps to the east. The southwestern part as well as the southeastern edge of the region show a local hail maximum with up to 3.1 hail days per year. The existence of these two hot spots may be explained by their proximity to the Mediterranean as during southerly flows, warm and moist air is advected preferably through the Rhone valley. The warm and moist air can then be lifted, for example, near a front system crossing the country from northwest to southeast, leading to forced convection. This effect was confirmed by Schemm et al. (2016), 365 who analyzed the relation between radar-based hail streaks over Switzerland and adjacent regions and cold fronts identified in high-resolution model data (COSMO-2;Steppeler et al., 2003, Jenkner et al., 2010  Midi-Pyrenees regions are the two regions well known in the literature for their high hail probability (Vinet, 2001, Punge et al., 2014. Hermida et al. (2015) used data from the ANELFA (Association Nationale d'Etude et de Lutte contre les Fléaux Atmosphériques) hailpad network and found that the Gers Department, located on the west side of the Midi-Pyrenees region, is the area the most affected by hail in southwestern France. The western and northern sides of the Pyrenees are also frequently 375 affected by hail with up to 2.5 hail days per year. According to Berthet et al. (2011), hail in that region frequently occurs when a low-pressure system is located over the western part of Spain leading to southwesterly flow over France associated with the advection of warm and moist air over the Pyrenean mountain range.
In Germany, the main hail hotspot is located in the Southwest in the federal State of Baden-Württemberg (region 4 in Figure   1), specifically over the Swabian Jura (mentioned as B on Figure 1) (Kunz and Puskeiler, 2010;Koebele, 2014). Moreover, Kunz and Puskeiler (2010) hypothesized that the southwesterly flow meets the Swabian Jura at a very sharp angle, which reduces the Froude number considerably and align the wind parallel 390 to the mountain chain. This flow modification is assumed to be responsible for the flow convergence at low levels as was also found in model simulations using COSMO-DE by Koebele (2014).
Another local maximum of up to 2.6 hail days per year is found North of the Alps, on the western part of the State of Bavaria (region 5 in Figure 1). This result is in good agreement with the conclusion of Nisi et al. (2018) who found that this region can be affected by around 3 hail days per year (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014).

395
In the northeast of Germany, a local maximum of up to 3.2 hail days per year is positioned over the Saxon Ore Mountains (labelled A in Figure 1) South of the city of Dresden. Note, however, that this maximum is mainly caused by a high number of SCS in the year of 2007 (Piper, 2017), which was characterized by frequent upper air troughs over Western Europe and ridges over Central Europe (Wernli et al., 2010), leading to high-pressure gradients on the eastern part of Germany in combination a southeast-to-northeast flow regime from the Czech Republic (note that the almost same situation occurred in 2019).

400
The northwestern part of Germany, including the States of Hesse (region 2 in Figure 1) and Rhineland-Palatinate (region 3 in

Annual variability
The frequency of SCS shows a very large annual and multi-annual variability (e.g., Nisi et al., 2018). This variability is partly related to large-scale flow mechanisms such as the presence of specific northern Hemisphere Teleconnection patterns representing the low-frequency mode of the climate system (e.g., North Atlantic Oscillation, NAO, or East Atlantic pattern, EA) or by variations in the sea surface temperature (Piper et al., 2019). Having reconstructed a very large event set of SCS/hailstorms 410 as presented in the previous section, we are also interested how the frequency of these events vary across the whole domain and regionally.
Averaged over the entire investigation area, the annual number of hail days is between 72 (2010) and 103 (2006) with a mean of 86 ( Figure 5). In 2006, large parts of Europe, including Germany, Belgium, Luxembourg, and northwest France, experienced higher temperatures than on average, especially during the end of June and July (NOAA, 2007), where two (moderate) heat 415 waves occurred (Fouillet et al., 2008). As a result, the sea surface temperature over the Mediterranean showed a positive anomaly (NOAA, 2007, Lenderink et al., 2009, leading to intense evaporation rates and, consequently, to an increase in the amount of water vapor in the atmosphere (Chaboureau et al., 1998). The spatial distribution of hail days in 2006 ( Figure 6) strongly resembles the climatology, with several maxima near hilly terrains and minima near the coastlines. Some hot spots can also be detected over the northwest part of France and in southwestern Germany.

Seasonal and diurnal development of SCS
The large spatiotemporal variability of hail discussed in the previous sections leads us to the question of the seasonal and diurnal development of SCS at the regional level. For this purpose, the entire study area is divided into 13 subdomains of similar 430 size (around 75,000 km 2 ) framed in Figure 3. We selected five subdomains with different terrain and climatological character-  To quantify the number of hail days in each subdomain, the average number of hail days for consecutive 10-day periods was calculated for the period 2005 to 2014 (Figure 7). Despite the large variability seen in the seasonal cycles of the subdomains considered, some similarities can be recognized. All time series of the different subdomains feature a clear annual cycle with a 440 minimum of hail days in spring and autumn and a maximum during the summer. This characteristic cycle with a strong increase in the hail day frequency during April/May, a significant decrease around September, and with a maximum during the summer months was found by several other authors such as Dessens (1986) and Vinet (2002) for France, Belgium and Luxembourg, Gudd (2003), Deepen (2006, Mohr and Kunz (2013) and Puskeiler et al. (2016) for Germany, and Nisi et al. (2014) andNisi et al. (2018) for Switzerland and northern Italy.

445
The mountainous subdomain MAS shows the largest average number of hail days and has the most pronounced annual cycle. Until the end of April, the average number of hail days for the 10-day running mean is below 10. During May and beginning of June, the number increases substantially from 9 around May 5 up to 17 days on June 4. The more pronounced diurnal temperature cycle for continental regions, associated with a higher lapse rate in combination with orographic lifting, may explain this increase (Berthet et al., 2011). After June 4, the average number of hail days increase steadily until to reach 450 the overall maximum for the MAS region at the end of July with 23 days.
Subdomain IDF likewise show a high hail frequency during the summer with up to 14 days, mainly at the end of July. This subdomain is under the influence of the Atlantic Ocean (Cantat, 2004), leading to an increased frequency of troughs (Vinet, 2001. Within this subdomain, the number of hail days increases slightly until the peak with a first local maximum in the beginning 455 of June (around 10 hail days) and a second local maximum at the beginning of July with around 12 hail days. Spring hailstorms may be associated with subtropical air masses coming from Spain, while summer storms preferably form ahead of cold fronts (Berthet et al., 2011). The number of hail days decreases sharply from the hail-peak season toward the end of September.
Subdomain SWF has a very broad hail peak in the middle of June with 12 hail days centered around June 14. Afterwards, the number slightly decreases and reaches 4 hail days at the end of September. This maximum found in June differs from the 460 analysis of Dessens et al. (2015), who found that May is the most active month followed by July over the southwestern part of France and the Mediterranean area (situated along the Rhone valley). Also Fraile et al. (2003) and Hermida et al. (2013) found that May is the month with the highest hail kinetic energy in southwestern France. Reasons for this discrepancy can be due to a longer period analyzed by Dessens et al. (2015), while Hermida et al. (2013) and Fraile et al. (2003) focused on a time range starting from the 90s.

465
Subdomain BAV, located in Southeast Germany, has the maximum of hail days at the end of July, later in the year than the other subdomains. Kunz and Puskeiler (2010) and Puskeiler (2013) also found that July is the month with the highest number of hail days in central and southern Germany. The development of hailstorms shown in Figure 8 represents the times where the CCTA2D detects the first radar reflectivity of 55 dBZ or more. Since the local time (LT) varies through Europe with approximately one hour from Brittany in France to Saxony in Germany, all times originally given in UTC are converted to LT, representing four minutes per degree of longitude.
In all subdomains, hail occurs most frequently in the afternoon between 13 and 18 LT, while between midnight and 10 LT the 475 fewest events are detected (Figure 8).
Some discrepancies appear in the daily cycle, mainly depending on the location and characteristics of the respective subdomain. For example, the frequency of hailstorms in BEL reveals a large increase during the afternoon (14-15 LT) and a slow, but gradual decrease toward the morning.
In contrast to the subdomains located in the northern part of Europe, domains MAS over the Massif Central and SWF slightly 480 peak one hour later at around 16 LT. The peak during the late afternoon for more continental regions is presumably due to local orographic effects, such as slope or valley winds (Nesbitt and Zipser, 2003).
For subdomain SWF, the average number of hail days remains high in the late evening (20 to 22 LT). 20 Figure 9. Histogram of all SCS mean length.
A plausible effect is that severe storms may develop from pre-existing scattered thunderstorms that form during the afternoon as was found by Nisi et al. (2016) and Nisi et al. (2018). This feature might be decisive for the hailstorm maximum in the 485 evening in the canto of Ticino in southern Switzerland.
Some literature exists regarding the diurnal cycle of hail in Europe (Punge and Kunz, 2016). Bedka (2011), for example, recognized a diurnal cycle of overshooting tops that is related to the presence of orography and/or to the distance to the sea. Kaltenböck et al. (2009) found a peak in hail occurrence in the middle of the afternoon through Europe. Kunz and Puskeiler (2010) identified for southwestern Germany a maximum in the number of damaging hail events during 13 and 18 LT. The 490 same peak is found in Alpine regions such as Italy (Morgan, 1973) or Switzerland (Nisi et al., 2016), where the maximum of hail occurs in the late afternoon and the minimum in the morning according to radar data analysis. Lukach et al. (2017) demonstrated for southeast Belgium that hail falls mostly during 15-16 UTC, which is in accordance to the daily cycle in subdomain BEL that includes Belgium. In this study we found a peak of hail around 16 LT in subdomain SWF, while Mallafré

Main characteristics of hail tracks
In the following section, we explore the main characteristics of our sample of radar-detected hail tracks. The length is defined as the distance in kilometers between the start and the end of a track determined by CCTA2D, i.e., the period where a threshold of 55 dBZ is reached or exceeded. The distribution of the lengths shown in the histogram in Figure 9 approximately follows an 500 exponential function with a maximum for the first class. In general, the mean length (with standard deviation) is 41.5 ± 36.4 km with a median of 29.5 km for the entire investigated area. The tracks reconstructed for Germany have a mean length of 39.1 ± 33 km and a median of 27 km whereas in France, the mean length is slightly larger with 43.9 ± 39.8 km and 32 km for the median. In total, 43% of all recorded storms over Western Europe have a length between 1 and 10 km. The number of tracks having a length between 10 and 20 km decreases to 19% of 505 the overall sample. Approximately 30% of all tracks have a length between 20 and 100 km, and less than 10% are greater than 100 km (not shown). Longer tracks can be expected for highly-organized convective systems, such as MCSs or supercells.
Only a few authors have analyzed hail tracks characteristics in West Europe, and only very few studies based their investigations over a sufficiently long period. Puskeiler (2013), for example, investigated hail tracks lengths using 3D radar data in Germany during 2005 and 2011 and found a mean length of 48 km with a strong decrease for longer streaks, inducing a high 510 standard deviation of 46.7 km and a median of 40 km. Dessens (1986) found a mean length of 80 km for a small sample of 30 hailstorms in southwestern France. Note that Dessens (1986)   The distribution of the hail track duration (Figure 10) is in accordance with the length, and also decreases almost exponentially with a peak at 30 minutes. As for the other physical characteristics, long-lived swaths are rare: only 2.4% of all cells persist over 5 hours (Fluck, 2018). The width, expressed as the maximum diameter of the largest reflectivity core (Z ≥ 55 dBZ) Figure 11. Same as Figure 10 but for the mean orientation. during a hailstorm, or the longest distance between two cores evolving laterally in cases of cell-merging or splitting, show approximately a Gaussian distribution. The peak is between 9 and 10 km (41% of all events; not shown).

520
However, as only the largest width of each swath is recorded, the results may be overestimated. The track angle shown in Figure 11 represents the mean orientation of a storm track, and is the angle recorded by CCTA2D at the center of a swath, as most of the storm tracks are approximately linear. The orientation is defined as the angle between the line intersecting the reflectivity core centers, before and after the central point of a swath, and the parallel of latitude intersecting the start point.
With this method, three time steps (i.e., 10 minutes) must exist for a track to be recorded. The maximum occurrence (19.3%) is 525 found in the propagation direction between 220 and 240 • , i.e., with a southwest direction. Around half (51%) of the hailstorms come from between 200 and 260 • . Only 2.8% of the swaths have a northwesterly direction. The principal southwesterly swath orientation found in the statistics confirmed previous results, such as that in the study of Puskeiler (2013), who found the highest number of hail days in southwest Germany with swaths oriented in southwesterly direction. In the Aquitaine region in France, Berthet et al. (2013) found that between 1952 and 1980 severe hailfalls came from a southwesterly direction with a 530 mean angle of 241 • .

Conclusions
This paper presents the first high resolution, radar-based hail statistics for a large central European region covering the countries of France, Germany, Belgium, and Luxembourg over a 10-year period. A European radar composite has been created from French and German national radar composites with a five minute time step and corrected with lightning data. Tracks of SCS 535 have been reconstructed using the tracking algorithm CCTA2D. Only grid points exceeding the Mason (1971) criterion for hail (Z ≥ 55 dBZ) have been used for the hail assessment. From the spatial analyses of the hail signals, the following main results are obtained: -The frequency of hailstorms shows a very large spatial variability across the investigation area. In general, there is a coast-to-continental increase in the number of hail days. While only a few radar-derived hail days occurred in the 540 northern parts of Europe (0-2 hail days per year), the number of hail days farther inland is much higher.
-Most of the hail hot spots are found on the leeward side of low-mountain ranges such as the Massif Central in France or the Black Forest in Germany. The large number of track onsets around orographic structures suggest a strong relationship between hailstorm occurrence and flow conditions induced or invigorated by orography such as a flow-around regime with subsequence flow convergence on the lee side.

545
-On the regional-scale, some differences in the seasonal and hourly distribution of hail occurrences are found across Europe: In Southwest France, for instance, the hail maximum is in mid-June, but occurs two months later in August in eastern Germany.
Our radar-derived hail frequency estimations and maps have, of course, several limitations and uncertainties. First, due to the local-scale nature of hailstorms and the lack of accurate observations, the reconstructed streaks and their statistics are 550 difficult to validate. No homogeneous monitoring system for hail exists over the entire investigation region, but only some local networks, for example, the hailpad network over the Southwest, Central, and Southern France operated by ANELFA are available. However, hailpad networks do not exist in Belgium, Luxembourg, and Germany. Thus, the use of radar and lightning data only provides proxys for hail occurrence.
Because only 2D radar data were available for this study, more sophisticated hail detection algorithm, such as those based on 555 echo-top height (e.g. POH) or vertical integrals of reflectivity (e.g. MESH), which generally show higher skill in hail prediction (e.g.Skripniková andŘezáčová, 2014, Kunz andKugel, 2015, Puskeiler et al., 2016), could not have been applied. The radar coverage over Western Europe is reliable, but several regions are still not or not sufficiently covered by several radars, such as the Alpine chain or some areas in the southeastern part of Germany and near Lake Constance. This leads to some data gaps in the final composite.

560
Despite the above mentioned limitations in our methods, the final results are in good accordance to other studies such as those for Germany based on 3D radar data (Puskeiler, 2013, Puskeiler et al., 2016, Schmidberger, 2018. The spatial distribution of hail signals in our study area is also similar to satellite-estimated hail frequency based on overshooting cloud tops as described by Punge et al. (2014) and Punge et al. (2017).
The investigations can be improved by extending the observation period until today. This is important especially in the sub-565 domains highly exposed to hail. More accurate detections of hail can be achieved via the use of dual-polarisation measurements (Heinselman and Ryzhkov, 2006). Furthermore, detailed investigations of the flow characteristics depending on atmospheric conditions, for example by using high-resolution numerical weather prediction models, can help to find robust evidence of the flow-around regime that may be decisive for the increased hail frequency downstream of several low mountain ranges, and can also contribute to a better understanding of the influence of orography on the triggering of convection.

570
Data availability. 2D radar data for Germany can be accessed via the German Weather Service (DWD) ftp server, while French national radar composites are available upon request to the French Meteorological Service (Météo-France). SCS/hail tracks were computed on radar data and are available upon request to M. Kunz. ERA5 data can be downloaded from the ECMWF server.
Author contributions. EF edited most parts of the paper, performed the statistical analyses and computed the SCS/hail tracks. MK verified in details the analytical methods and results, added crucial suggestions to the paper and contributed to the editing/revision of the manuscript.

575
PG and SR added constructive suggestions to the paper. MK supervised the project in collaboration with PG and SR.
Competing interests. The authors declare that they have no conflict of interest.