Daytime identification of summer hailstorm cells from MSG data

Identifying deep convection is of paramount importance, as it may be associated with extreme weather phenomena that have significant impact on the environment, property and populations. A new method, the hail detection tool (HDT), is described for identifying hail-bearing storms using multispectral Meteosat Second Generation (MSG) data. HDT was conceived as a two-phase method, in which the first step is the convective mask (CM) algorithm devised for detection of deep convection, and the second a hail mask algorithm (HM) for the identification of hail-bearing clouds among cumulonimbus systems detected by CM. Both CM and HM are based on logistic regression models trained with multispectral MSG data sets comprised of summer convective events in the middle Ebro Valley (Spain) between 2006 and 2010, and detected by the RGB (red-green-blue) visualization technique (CM) or C-band weather radar system of the University of León. By means of the logistic regression approach, the probability of identifying a cumulonimbus event with CM or a hail event with HM are computed by exploiting a proper selection of MSG wavelengths or their combination. A number of cloud physical properties (liquid water path, optical thickness and effective cloud drop radius) were used to physically interpret results of statistical models from a meteorological perspective, using a method based on these “ingredients”. Finally, the HDT was applied to a new validation sample consisting of events during summer 2011. The overall probability of detection was 76.9 % and the false alarm ratio 16.7 %.


Introduction
Measurements of solar reflection and emittance of cloud systems by means of satellite sensors have been shown to be in-strumental for retrieving the cloud optical and microphysical properties for a variety of uses, such as cloud physics, meteorology and climate studies (e.g., King et al., 1992).The EUropean organisation for the exploitation of METeorological SATellites (EUMETSAT) has established the Satellite Application Facility on Support to Nowcasting and Very Short Range Forecasting (SAF-NWC), which makes available the algorithms for retrieving cloud physical properties.Marcos and Rodriguez (2013) developed an algorithm that provides an estimate of the probability of precipitation occurrence using information on the microphysical properties of the cloud top.The radiative properties of a cloud were characterized by means of the effective radius (R e ) and cloud optical thickness (OT; see acronym list in Appendix A).
The current generation of European geosynchronous satellites yields a high-quality signal and enhanced spatiotemporal resolution, which represent a major step forward for monitoring of short-lived weather phenomena such as rapidlydeveloping convective storms, for which high spatial and temporal resolution is critical.The first such sensor is the Spinning Enhanced Visible and InfraRed Imager (SEVIRI) (Schmetz et al., 2002), which is the main instrument on board the European geostationary satellite Meteosat Second Generation (MSG).It has 12 spectral channels with spatial sampling distance of 3 km at the sub-satellite point and a high resolution visible (HRV) channel with spatial sampling distance of 1 km.Temporal resolution for the full disk of the SE-VIRI is 15 min, with the possibility of obtaining rapid scans at shorter time intervals.
The increased number of spectral channels of SEVIRI over the previous generation Meteosat sensors makes it possible to develop different multispectral and multithreshold techniques to identify cloud types, such as the split window technique (Inoue, 1987).For example, when the BT A. Merino et al.: Identification of summer hailstorm (brightness temperature) difference of channels at 11 and 12 µm is greater than 2.5 K, the cloud is considered cirrus (Kurino, 1997).Other authors (Strabala et al., 1994) use BTs in the spectral range of 8-12 µm to identify the cloud thermodynamic phase.This trispectral technique is based on the fact that the absorption coefficient for water particles increases more between 11 and 12 µm than between 8 and 11 µm; for ice, the reverse is true.
Many studies have focused on identification of storm cells using various satellite data.Kurino (1997) found that the BT difference of 11-6.7 µm is 0 K or less for convective clouds associated with heavy rain.Schmetz et al. (1997) found that the equivalent BT of the 6.7 µm channel can be larger than that of the 11 µm channel by 6-8 K.This is because deep convective clouds penetrate the stratosphere, injecting water vapor there.The temperature in the stratosphere is warmer than that in the upper troposphere, so it is often true that BT in the water vapor channels is higher than BT in the thermal infrared channels.Zinner et al. (2008) used the temperature index in the tropopause, obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) model, to detect mature convective clouds.The authors found that these clouds had 1.5 K TB at 6.7 µm less than the tropopause temperature.Cattani et al. (2009) studied the cloud optical and microphysical characteristics of convective storms.They found that clouds with precipitation intensities greater than 5 mm h −1 had an effective radius (R e ) between 20 and 30 µm at their tops.Other authors have studied microphysical characteristics of convective storms (Rosenfeld et al., 2008;Mecikalski et al., 2011).Updraft speeds could be computed at cloud tops during the developing phase using satellite sensors, according to the cloud-top cooling rate given a high temporal resolution.Using the "rapid scan" mode of geostationary satellites it is possible to analyze microphysical properties of convective cloud growth (Mecikalski et al., 2011).Additionally, Rosenfeld et al. (2008) conceived a method to infer the intensity of the storm updrafts from the microphysics of cloud-top particles.Vertical profiles of R e were computed from the study of areas with convection in various stages of development.The authors concluded that cloud tops with very strong updrafts contain small ice crystals, and reflectance in the near-infrared (NIR) channels is high.
Another feature of severe storms (associated to severe weather according to Johns and Doswell, 1992) is that they often develop overshooting tops with a V-shape leeward of the cloud top, resembling a diverging plume above the anvil top (Heymsfield et al., 1983).This plume can have high reflectance in the NIR channels because it contains small ice particles (Levizzani and Setvák, 1996).In fact, Adler et al. (1983) found that most storms with this V-shape were related to severe weather (tornadoes, hail, and intense rain).The presence of overshooting clouds confirms strong updrafts within the storm often being associated with hazardous weather (Dworak et al., 2012).These types of clouds can be detected using spatial gradients of 11 µm BTs (Bedka et al., 2010;Bedka, 2011).Setvák et al. (2010) observed with enhanced infrared window satellite imagery that deep convective storms can have long-lived cold rings at the cloud top, causing a warm area inside the ring in overshooting clouds.
In summary, numerous studies have identified different aspects of convection using satellite data.The strong relationships established between different products and hail precipitations have been considered in preparing this study, introducing an unsupervised objective hail detection algorithm to identify hail precipitation in real time.
The middle Ebro Valley (MEV) in the northeast of the Iberian Peninsula is one of the areas in Europe with highest frequency of hail events, with about 60 days characterized by storms each summer (López and Sánchez, 2009).Since 2001, the Atmospheric Physics Group (GFA) of the University of León in Spain has been developing a number of projects in this area to study hailstorm convection and monitor its development (García-Ortega et al., 2007, 2011;Sánchez et al., 2009).The GFA uses a C-band weather radar system with a nowcasting model for detection of hailstorms (NMDH) (López and Sánchez, 2009).However, this hail detection system has some drawbacks, such as limited spatial range and radar beam shielding in mountainous areas.The aim of the present study is to develop a tool for identifying hail storms in real time that avoids the drawbacks of the radar system.Thus, the aim was to develop a nowcasting tool to identify hail-bearing clouds using MSG data.Its temporal resolution (15 min in operational mode) is lower than that of the radar (4 min for the University of León radar), but it may be used to monitor convection in real time.
To meet this objective, a hail detection tool (HDT) has been developed in two steps using logistic regression models.First, the deep convection is identified using a convective mask algorithm (CM); second, the hail mask algorithm (HM) is used to identify hail precipitation within the clouds.This system was applied in the summer months (June, July and August) during daylight hours (with solar zenith angle lower than 70 • ), when hailstorms are most frequent in the study area.

Study area and nowcasting model for detection of hailstorms
The study area covers the entire Iberian Peninsula (Fig. 1), although the ground truth data (hail precipitation intensity) were obtained exclusively for the MEV.As already noted, the GFA uses a C-band weather radar system with an NMDH (López and Sánchez, 2009).The NMDH was based on a database of 702 events, identified by the radar using Thunderstorm Identification, Tracking Analysis and Nowcasting (TITAN) software.The NMDH was trained and validated using a network of 700 voluntary observers across most of the MEV, providing a probability of detection (POD) rate of 85 % and false alarm ratio (FAR) of 15 %.Therefore, this tool determines the presence or absence of hail and predicts the spatial likelihood of hail precipitation for each storm detected by the radar.
In this paper, the results of the NMDH have been used as "ground truth" for the construction of the training and validation databases of HDT.The reason for this is the accuracy of this tool in identifying hail precipitation (López and Sánchez, 2009).Furthermore, it is possible to determine to some extent the time of hail precipitation events.Given the small spatial and temporal scales of these events, determining the exact time is important.Nevertheless, it was decided not to use direct observation of hail data on the ground from the observer network, as these reports have larger time uncertainties.

The logistic regression method
Logistic regression is a widely used statistical tool in meteorology and land use studies (Applequist et al., 2002;Bastarrika et al., 2011;López et al., 2007).Sections 4.2 and 5.2 describe construction of two logistic regression models for formulating the CM and HM algorithms.The present section describes basic features of the logistic regression technique applied in the study.The technique provides probability of occurrence (P) of a particular weather event (categorical variable) from values of a number of metric explanatory variables (X k ).When the categorical variable is dichotomous, binary logistic regression is used.When there are more than two values involved, the latter is substituted by multinomial logistic regression (Hosmer and Lemeshow, 1989).The P obtained by the model for a categorical variable of the two groups is expressed as follows: where β j X j ; X j for j = 1, . .., k are the metric explanatory variables; k is the total number of variables or interactions between variables included in the model; α is the intercept; β j for j = 1, . .., k the various discriminating coefficients; and P n represents the probability of belonging to group (1) or group (2), which ranges between 0 and 1.Note that p 1 + p 2 = 1.
As previously indicated, logistic regression was used to generate the CM and the HM.In both cases, the metric explanatory variables are the MSG channels, except the HRV (Table 1).For construction of the CM, the categorical variable was the presence or absence of cumulonimbus and, for the HM, the presence or absence of hail.The model is constructed using the stepwise method, where the metric explanatory variables can be added according to their forecasting potential (Hair et al., 1999).There are different methods for incorporation of the variables into the model: a higher Wald coefficient, higher conditional probability, or the likelihood ratio.In this case, the criterion of the likelihood ratio was used, which is the greatest reduction in the value of −2 log likelihood (−2LL).Here, the variables were introduced when the reduction of −2LL had significance level 0.05 and, for the output, significance level 0.1.The predictive equations were constructed with purely statistical criteria.However, the combination of variables must have a meteorological justification for the analyzed weather phenomenon (Doswell and Schultz, 2006).According to Doswell et al. (1996) the "ingredients-based methodology" is the scientific path to assess an equation formulated using a strictly statistical method.
The logistic regression model can be improved by introducing interactions between the predictor variables, thereby increasing the predictive power of the equations.In this study, the model was constructed by introducing first-order interactions between the metric explanatory variables.Thus, data from a particular MSG channel may be input individually or combined with those from another channel.The presence of interactions causes the relation between the categorical variable and metric explanatory variable to depend on the value of a third variable.Thus, the contribution of each variable in the logistic model is conditioned by the variable with which it interacts.This fact complicates interpretation of the variables included in the model; however, it increases its predictive power significantly.The sign of β coefficients determines the sign of the contribution of each variable to the model, when there are no interactions, or also the interaction variable as a whole (X m ).In this case, a positive β coefficient involves greater probability in the model when the value of the variable increases, and the reverse is also true.When the variables show interactions they must be interpreted together physically, because in order to extract the sign of each single variable (X j ) it is necessary to combine their β coefficients (Eq.3), which depend on the value of the variable with which it interacts (X i ): where X j and X i are single variables involved in the interaction, and Sign X j is the sign of the contribution of this variable to the model.X m is the interaction variable (X m = X j X i ) and β m is the interaction coefficient.As the interaction variable (X m ) has a high order of magnitude, it is necessary to provide the interaction β coefficient with a sufficient number of decimal points.
To check whether the explanatory variables added to the model have a high degree of statistical significance, global fitting of the model may be computed via the chi-squared test for variation of the value −2LL with respect to the base model (without variables).Thus, a model with good fit will have a large reduction for −2LL, and a perfect fit is one in which the likelihood is 1 and −2LL is zero.Owing to the introduction of the independent variables, the chi-squared test for assessing significance of the reduction of −2LL must be significant.In addition, several different R 2 measures have been constructed to represent a global fit of the model.Some of these measures are the Cox and Snell (Hair et al., 1999) and Nagelkerke R 2 values (Hair et al., 1999); both are between 0 and 1, and 1 is considered a perfect fit.
Finally, the Wald test was applied to assess the significance of the β coefficients for variables included in the logistic equation (Hair et al., 1999).Variables with significance of 0.1 may thus be interpreted as metric explanatory variables in discrimination of the categorical variable.It is important to point out that the greater the significance and Wald statistic are, the greater weight the variable will have in the logistic model.The Wald statistic is used to compute the weight of each variable in the model instead of the β coefficients, because the metric explanatory variables have different orders of magnitude.
In accord with Hosmer and Lemeshow (1989), the cutoff point for discrimination between Eqs. 1 and 2 was fixed at 0.5, for both the CM and HM.

Convective mask
HDT consists of two distinct stages.First, the CM was set up to identify deep convection using MSG data.To do this, a database of cloud events was created and subsequently analyzed from a microphysical point of view.These microphysical and optical analyses of the various cloud types were studied to interpret the results of the logistic model from a meteorological perspective.Although there are numerous CMs for SEVIRI (Berendes et al., 2008;Henken et al., 2011), we opted to develop a new CM to verify the HDT in its entirety.However, other CMs for SEVIRI can be used as a basis for later application of the HM.

Database
The training database used to construct the CM included satellite imagery from the summer months between 2006 and 2010 during daytime (solar zenith angle lower than 70 • ).Cloud types were identified using red-green-blue (RGB) combinations.These combinations enabled representing physical information in the MSG channels (Lensky and Rosenfeld, 2008).Using the RGB compositions "day microphysical", "day solar" and "convective storms" by Lensky and Rosenfeld (2008), 700 events were classified, 100 for the following cloud types: cumulonimbus, stratocumulus, cirrus, nimbostratus, cirrostratus and multilayer cloud, stratus, and clear sky.An event was defined as an MSG image when cloud systems of the above-mentioned cloud types with a spatial extent of at least 10 km 2 could be identified by RGB combination.Only independent events were taken into account, removing MSG images related to the temporal evolution of the same cloud system.
The combination of physical properties of the clouds that make up these RGB schemes (cloud drop size, thermodynamic phase, cloud top height, optical thickness) permit each cloud type to be identified, following guidelines established by Lensky and Rosenfeld (2008)."Day microphysical" compositions (red = 0.8 µm; green = 3.9 µm reflectance; blue = 10.8 µm) are used to observe the microstructure of water and mixed-phase clouds during daylight.In this scheme, stratus with small cloud drops and high optical thickness appear white, whereas warm precipitating clouds with large particles are pink.Cirrus clouds with tops composed of large ice particles appear dark.The "day solar" scheme (red = 0.8 µm; green = 1.6 µm; blue = 3.9 µm reflectance) is very sensitive to cloud microphysics and serves to determine the particle size and phase at the cloud top.Finally, the "convective storms" scheme (red = BT difference 6.2-7.3 µm; green = BT difference 3.9-10.8µm; blue = reflectance difference 1.6-0.6 µm) highlights clouds with very cold tops and allows better identification of young, severe storms.Thus, convective clouds characterized by strong updrafts that appear bright yellow can be distinguished from cumulonimbus clouds with large ice particles, which show up as red (Lensky and Rosenfeld, 2008).Uncertainty in the detection of each type of cloud via these schemes has been reduced by combining information from the various RGB schemes used.
For the total 700 identified events, radiances and reflectances of the MSG channels were retrieved for each event at 15 min temporal resolution .The reflectances were transformed into albedos following the Lambert method to avoid their dependence on solar zenith angle, and the radiances were transformed into BTs.The 3.9 µm radiances were converted to BT, taking both contributions into account.It is important to point out that surfaces with high reflectances in the direction of the satellite may have albedos greater than 100 %.Only albedos and BTs of different MSG channels were then used as inputs in building the logistic regression model.

Results of logistic regression
As mentioned in Sect.2, logistic regression was used to construct the CM, with data from the first 11 MSG channels chosen as explanatory variables (Table 1) and the categor- ical variable represented by the presence (P 1 ) or absence (P 2 ) of cumulonimbus clouds.Thus, albedos and BTs of 600 cumulonimbus-free events were input to the model, along with the cumulonimbus events, to find the combination of MSG channels that best distinguished these two groups statistically.The binary logistic model is based on the forward stepwise method and executed over 10 iterations, introducing 9 variables in Eqs. ( 1) and ( 2).The estimated β coefficients gave a correlation coefficient equal to zero, according to Wald's test (Table 2); thus, all are susceptible to physical interpretation.Global fit of the model was assessed using the statistics described in Sect. 2. The chi-squared test used to assess the reduction of the −2LL parameter showed that it was significant (Sig < 0.05), so the fit was satisfactory.In addition, the Nagelkerke's R 2 and the Cox and Snell's R 2 gave values of 0.456 and 0.804, respectively, which can be considered an acceptable fit (Hair et al., 1999) (Table 3).Table 4 shows the classification of events included in construction of the logistic equation.A contingency table was built, and 700 events were classified according to the algorithm result.Of the 600 noncumulonimbus events included in the training set, only 9 were classified as cumulonimbus.Of the 100 cumulonimbus events, 12 were wrongly classified.The correct classifications were 97 %, with POD of 88 % and a FAR of 9 %.

Physical interpretation of results
To assist in physical interpretation of the logistic model results, the cloud physical properties optical thickness (OT), effective radius (R e ) and liquid water path (LWP) were extracted via the visible infrared solar-infrared split-window technique (VISST; Minnis et al., 1995).Data were obtained from the NASA Langley Research Center Cloud and Radiation Group website (http://angler.larc.nasa.gov/)with the same temporal and spatial resolution of MSG data.The microphysical properties dictate the emittance and bidirectional reflectance of clouds, since they determine the amount of radiation transferred to the surface, reflected and scattered to the satellite sensor, and the amount of radiation absorbed by the cloud layer.These characteristics vary according to cloud type (King et al., 1992).

Cloud physical properties: liquid water path, optical thickness and effective radius
The cloud physical properties were examined for a representative case study, for example, from the event on 12 August 2011 at 14:00 UTC (Coordinated Universal Time).This is because at that time there were a variety of cloud types over the Iberian Peninsula.Figure 2 is the RGB image "day solar", showing areas with stratus, nimbostratus, cirrus and cirrostratus clouds.To the right of the figure, the RGB image "convective storm" highlights cumulonimbus clouds in different stages of development.
The cloud physical properties were obtained for each cloud pattern and are described below.

a. LWP
The scatter plot of LWP values as a function of BT (10.8 µm) for the various cloud types (Fig. 3, upper panels) shows that cirrus and cirrostratus have the smallest LWP.Cumulus clouds with little development also had low LWP values, whereas stratus and nimbostratus are characterized by LWP values of up to 410 g m −2 , despite their relatively small thickness.Convective clouds showed the highest LWP, with many pixels having more than 1000 g m −2 .

b. OT
The OT (Fig. 3, middle panels) shows that cirrus clouds had very low OT values, since they are nearly transparent.OT values increased for cumulus clouds with weak development and, in stratus and nimbostratus clouds, they neared 50, since these clouds are opaque.Both developing and dissipating cumulonimbus clouds had the highest OT values (128 in this example), with no major differences between the two.

c. R e
The R e scatter plot (Fig. 3 bottom panels) shows that cloud tops below the glaciated zone had low R e values, as with stratus and developing cumulus clouds.In contrast, cloud tops containing ice crystals had much higher values.For strong convection in the developing stage, R e values were not particularly high compared to convection in the dissipating stage.This is due to the fact that updrafts within convective clouds in the developing stage are much stronger and, as a result, heterogeneous nucleation produced in the mixed phase does not have enough time to form large ice crystals.Thus, smaller ice particles are formed via homogeneous ice nucleation above the level of glaciation compared to particles of greater size formed via heterogeneous nucleation (Rosenfeld et al., 2008).
The microphysical properties of the different cloud types were investigated to determine features that enable distinction between convective and other cloud events.Cumulonimbus clouds are characterized by high LWP, high OT values and, depending on their stage of development, variable R e values.Thus, albedos and BTs included in the logistic model should reflect these characteristics.

Model input variables
The SEVIRI spectral channels and their combinations that were included in the logistic equation are shown in Table 2. Later, a physical interpretation of the variables, with respect to their contributions to the model, was developed according to the Wald statistics (Table 2).Variables showing interactions were interpreted together as follows.-Channels at 6.2 and 7.3 µm, and their interaction (channel 6.2 µm × channel 7.3 µm): emission in this part of the spectrum occurs within the water vapor (WV) absorption band.The channel at 6.2 µm is sensitive to WV emittance in the upper troposphere, and that at 7.3 µm to emittance in the mid-troposphere.These channels are also important when computing atmospheric corrections (Minnis et al., 1995).The 7.3 µm channel and its interaction with the 6.2 µm channel are the greatest weight in the equation (Wald parameter in Table 2).Cumulonimbus clouds are the only ones often extending through the entire troposphere, and this is a crucial feature for identifying this cloud type.These clouds produce high LWP values and low BTs in the WV channel spectrum.Together, these three variables are inversely correlated with high probabilities of cumulonimbus, so low BTs in the 6.2 and 7.3 µm channels are associated with high probability of cumulonimbus, which accords with physical expectations.
-Channels at 8.7 and 1.6 µm, and their interaction (channel 8.7 µm × channel 1.6 µm): the channel with the next-highest weight in the model is of 8.7 µm (Wald parameter in Table 2).The interpretation of this channel in the model has been carried out considering its category (infrared windows).This means that introducing new channels of the same category does not lead to a statistically significant improvement in the model.This fact was checked by substituting the 8.7 µm channel with other infrared window channels (at 10.8 and 12 µm).The results of this new model were very similar to the original, after adjusting the coefficient.The atmospheric window channels permit distinction between different cloud top temperatures.Developing cumulonimbus clouds have tops formed by ice crystals near the tropopause, and the channel therefore distinguishes them from middle and low clouds.The channel with which it interacts (1.6 µm), albeit with less weight in the model, permits distinguishing the cloud-top phase.At 1.6 µm, ice and liquid-water clouds have very different reflectances (Cattani et al., 2007;Rosenfeld et al., 2008), owing to sensitivity to R e and cloud phase.Cloud tops formed by liquid-water hydrometeors have lower absorption than those formed by ice crystals, so reflectances are higher for water than for ice particles.In this case, the signs of the variables are related.To extract these signs, it is necessary to combine their β coefficients (Table 2) using Eq. 3. Thus, channel 1.6 µm has a positive effect for channel 8.7 µm BTs of less than 226.91 K, and a negative effect for BTs greater than that value.This result has a physical interpretation, since clouds with 8.7 µm BTs of less than 226.91 K (lower than the level of homogeneous nucleation, according to Rosenfeld et al. (2008)) are formed by ice crystals.Greater albedo values at 1.6 µm mean that particles at cloud top are small, and updrafts are more vigorous (Rosenfeld et al., 2008); that is, the probability of cumulonimbus increases.On the contrary, clouds with 8.7 µm BTs warmer than this temperature may be formed by liquid water, with large albedos at 1.6 µm, which diminishes the probability of cumulonimbus.
-Channels at 3.9 and 0.8 µm, and their interaction ( channel 3.9 µm × channel 0.8 µm): the channel with the next-highest weight in the model is at 3.9 µm (Wald parameter in Table 2).BTs in this channel, together with those measured in channel 1.6 µm, are very sensitive to R e at cloud tops (Cattani et al., 2007;Rosenfeld et al., 2004).Cloud tops formed by water or small ice hydrometeors have low absorption and thus high BT.This fact enables distinguishing between clouds with tops characterized by low R e values and high TB in this channel, and those with tops associated with high R e values and low BT.This channel interacts with the 0.8 µm channel and, according to Nakajima and King (1990), reflectances at this wavelength depend primarily on cloud OT.As a result, this channel distinguishes between very dense, opaque clouds and thinner, transparent clouds.The interaction between these two channels couples the effects of R e and OT, and their information can be used for simultaneous retrieval of OT and R e (King et al., 1992).Mecikalski et al. (2010) showed that changes in the reflectances of these channels were related to the variations of OT and hydrometeor size in growing clouds.The sign of the variables in the model also varies.Combining their β coefficients (Eq.3), the 0.8 µm channel has a positive effect over the entire range of possible values of the 3.9 µm channel.In other words, a large albedo at 0.8 µm produces large OT and a greater probability of cumulonimbus.The 3.9 µm channel has a positive effect for 0.8 µm channel albedos greater than 125.44 %, and a negative effect for lesser albedos.The physical reason is similar to that explained in the second point above.Albedos greater than 125 % are typical of very dense clouds and, since BTs increase in the 3.9 µm channel, the probability of cumulonimbus increases.According to Mecikalski et al. (2010), reflectances in the visible channels do not provide definitive information about cloud top character, but they can be used to confirm that the cloud is optically thick and the use of 3.9 µm is acceptable.
As mentioned above, the CM uses a large number of variables within the discriminating equation.This is because of the wide range of cloud types included in the database.However, all channels selected using statistical criteria have an adequate physical interpretation in cloud-type discrimination.The results reveal that the "ingredients" required to distinguish cumulonimbus clouds from other cloud types are, in order of importance, as follows: -water vapor in the troposphere (assessed by the WV channels), -thermodynamic phase of particles at cloud top and top cloud temperature (assessed by the NIR channels and the 8.7 µm channel), and -cloud optical thickness and size of the particles at cloud top (assessed by the NIR channel and the 0.8 µm channel).
Cumulonimbus clouds are thick, with a high WV concentration throughout the troposphere, and their cloud tops are formed by ice crystals.The combination of all these features enables discriminating this cloud type from others, and so these are the three "ingredients" that must be included in an algorithm for detecting cumulonimbus clouds.These ingredients are linked to the aforementioned physical properties.WV concentration, for example, is related to LWP.The thermodynamic phase of the cloud is related to R e , since this parameter strongly depends on the presence of ice or water.

Hail mask algorithm
The second step in construction of the HDT from MSG data is identification of hydrometeors within cumulonimbus clouds detected by the CM.The method for building the HM is similar to that described above for the CM.First, a database of cumulonimbus events was created, with and without hail.
Then, the logistic regression model was constructed using the albedos and BTs to determine the MSG channels that best distinguish the hail events.Finally, microphysical and optical cloud properties were studied in areas with hail to interpret the MSG channels included in the model from a physical perspective.

Database construction
The database for the HM was built using information from the GFA weather radar.This database includes daytime events recorded in the summer months, during five observation campaigns between 2006 and 2010.To distinguish hailbearing from hail-free events, the results of the NMDH implemented for the radar data (López and Sánchez, 2009) were used.The NMDH provides the likelihood of hail precipitation, and its results were considered "ground truth" data.A hail pixel is one at which the likelihood of hail according to the NMDH exceeds 90 %.To build the database using radar data, the following issues were taken into account.
-Temporal resolution.The radar data must correspond to the same time span of the MSG satellite scan across the study area.In this case, for the MEV the satellite gathers data between 8.5 and 9 min after the scan initiation.
-Parallax effect.The parallax effect is important for the MSG because the satellite is above the Equator, taking measurements of Europe at relatively low angles.Lábó et al. (2007) found a deviation up to four pixels towards the southwest for high clouds over Hungary.The MEV is not situated at the satellite nadir and the analyzed cumulonimbus clouds have high tops; thus, correction must be done.The deviations computed in this study follow the method of Vicente et al. (2002), which takes data of cloud-top height from vertical reflectivity profiles obtained by the radar.For example, for high clouds (between 15 and 18 km) the deviations reached 18 km to the south; there were very small deviations in the E-W direction, because the zero degree meridian crosses the study area.Eventually, this correction enables comparing satellite and radar data at the same surface location.
-Spatial resolution.The GFA radar data have a resolution of 0.75 km, whereas the MSG data have spatial sampling distance of 3 km at the subsatellite point.Because of this and considering small deviations of hail precipitations that may be attributed to wind, we considered only hail precipitation with an extent of at least 18 radar pixels.In this case, an event is defined as a cumulonimbus with a spatial extent of at least 10 km 2 identified in radar images pertaining to independent convective cells.
The hail training database was thus constructed, with the following events: -100 events of precipitating cumulonimbus clouds with hail, -50 events of precipitating cumulonimbus clouds without hail, and -50 anvil clouds.
Finally, radiances or reflectances of the first 11 MSG channels were considered for each event (Table 1).These were transformed into BTs and albedos, respectively, as with the CM algorithm.

Results of logistic regression
The HM was built with a categorical variable of binary type, taking value 1 for hail events and 0 for hail-free events.The binary logistic model was created as with the CM, introducing the albedos and BTs of both hail-free and hail events.The model executed eight iterations, inputting four variables in the equation.The β coefficients estimated for each variable are significant (Sig < 0.1) for physical interpretation, according to Wald's test (Table 5).Global fit of the model was assessed through the statistics described in Sect. 2. The chi-squared test used to assess reduction of the −2LL parameter yielded significant results (Sig < 0.05), revealing a good fit for the model.In addition, Nagelkerke's R 2 and Cox and Snell's R 2 indicators yield 0.784 and 0.588, respectively, demonstrating an acceptable fit (Hair et al., 1999) (Table 6).Table 7 shows the classification of events included in the construction of the logistic equation.A contingency table was built, classifying the 200 events according to the model.Among the 100 hail-free events, only seven were wrongly classified.Among the 100 hail events, eight were wrongly classified.Correct classification amounted to 92.5 %, with POD 92 % and FAR 7 %.

Physical interpretation of results
As in the case of the CM, the cloud physical properties OT, LWP and R e were extracted from the VISST algorithm for a sample case study, to make a proper physical interpretation of the MSG channels selected by the model.The aim was to determine microphysical characteristics in areas with and without hail.To study these properties, several singlecell storms were chosen.The NMDH was used to identify storm areas with high likelihood of hail precipitation.Then, scatter plots of the cloud properties as a function of corresponding cloud-top temperatures were used to compare these hail sectors with others within the cumulonimbus cloud.

Cloud physical properties: liquid water path, optical thickness and effective radius
The scatter plots were constructed by selecting convective areas identified by the CM in the sample case study (14:00 UTC, 12 August 2011).Areas affected by hail are shown by black contours (Fig. 4).The radar only covers the northeast of the study area (red circle in figure), so hail areas outside of this coverage are not identified.
a. LWP The LWP (Fig. 4, upper panels) did not reveal major differences between hail and hail-free areas.Nevertheless, whereas clouds without hail and lower heights exhibit a wide range of LWP values of up to 3000 g m −2 , clouds with hail have smaller LWP values.For hail areas, most of the water vapor is transported by strong updrafts toward upper levels in the cloud, forming ice.However, in other parts the cloud is affected by downdrafts and, thus, there are significant accumulations of liquid water at the base, thereby increasing the LWP.

b. OT
The OT (Fig. 4, middle panels) shows that the hail areas were mostly in regions with OT between 40 and 100.It is seen that many cloud pixels not associated with hail precipitation also have high OT values, because the anvil has high values.
c. R e The R e (Fig. 4, bottom panels) shows that most hail pixels had R e between 30 and 50 µm, considerably smaller than values at other cold pixels.This represents the microphysical property that best discriminates hail and hail-free regions within a cloud.In hail areas, updrafts are stronger than in other parts of the cloud, and lower R e values are found.

Variables input to the model
Channels included in the logistic equation are shown in Table 5.A physical interpretation of the variables, along with their interactions in order of weight in the model according to the Wald statistics, is shown below.
-Channel at 0.8 µm: this channel with a positive effect has the highest weight in the model (Wald parameter, Table 5).Thus, the greater the 0.8 µm albedo, the greater the probability of hail.According to Heymsfield et al. (1983), hail areas are commonly associated with high OT, related to overshooting clouds or a V-shaped form.Bedka (2011) found that 53 % of cumulonimbus with overshooting clouds produce hail on the ground.To detect this type of structure it is necessary to apply techniques of spatial recognition, since the present method using a pixel-by-pixel analysis does not permit such detection.However, it has been observed that these structures have high albedos in the visible and NIR channels (Berendes et al., 2008).Thus, pixels with elevated albedos in the 0.8 µm channel must be considered for hail identification.How-ever, apart from the visible channels, additional channel information is necessary for this detection (Berendes et al., 2008).
-Channels at 6.2 and 1.6 µm, and their interaction (channel 6.2 µm × channel 1.6 µm): the channel with the second-highest weight in the model is WV 6.2 µm (Wald parameter, Table 5), along with its interaction with NIR 1.6 µm.Channel 6.2 µm has a negative effect for 1.6 µm albedo values of less than 56.96 % (this threshold is obtained by combining the β coefficients Thus in practice, channel 6.2 µm will have a negative sign and channel 1.6 µm a positive sign.The probability of hail increases with albedo in the 1.6 µm channel and BT diminishes in the 6.2 µm channel, consistent with physical expectations.Storms can produce hail in their developing stage when their mass centroids are in the upper layers (López and Sánchez, 2009).This results in high WV concentrations in the upper tropospheric layers, producing a decrease in BTs at 6.2 µm.Channel 1.6 µm does not have the greatest weight in the equation; however, it is fundamental for hail detection.Its reflectances are associated with ice particle size, since cloud tops formed by liquid water are filtered by the CM.
As seen in the analysis of microphysical cloud properties, hail areas have small ice particles at their cloud tops because of the strong updrafts.It is well known that presence of hail on the ground is directly related to updraft speed (López et al., 2000).This speed determines R e at cloud tops (Rosenfeld et al., 2008).Therefore, reflectances in this channel are higher than those in regions with large ice particles.These results reveal that the ingredients necessary for discriminating regions with hail precipitation within cumulonimbus are, in order of importance, as follows: -optical thickness (assessed by the visible channels), -water vapor in upper troposphere (assessed by the WV channel), and -speed of updrafts (assessed by the NIR channels).
The combination of all these ingredients discriminates hail sectors from the remainder of the cumulonimbus cloud.

Verification of HDT from MSG data
The HDT was verified with convective hail events recorded in the MEV during summer 2011, since these data were not included in the initial training database.

Database
The same procedure as that chosen to extract hail events using the NMDH (no direct ground measurements) was followed to build the verification database.Hail-free events were either associated with convective clouds with tops higher than 10 km and with liquid precipitation registered on the ground of various intensities (radar-measured), or with parts of the anvil of the convective cloud.The verification database includes the following events: -26 convective-free events (cirrus, cirrostratus, stratocumulus, stratus and blue-sky), -26 hail events corresponding to convective cells (one or more cells), and -26 hail-free events (rain of varying intensities and anvil).
The database constructed using these radar data was considered ground truth.To assess the HDT results, we need to compare probabilities obtained by the model built using MSG data and the ground truth data for each event.Spatial probability weighting was used for this comparison.For each event, the central MSG pixel was extracted together with eight surrounding pixels, and the maximum likelihood among all nine MSG pixels was considered.The reason for this spatial weighting is that data from the MSG pixel might not coincide exactly with hail recorded on the ground.Deviations may occur, one source of these being the error in cloud-top estimation with the radar, and another the computation of the parallax effect.Other deviations may be from strong wind shear, which can tilt the storm.In all these cases, the area on the ground where hail is recorded does not coincide exactly with the cloud top.

Results
The results were assessed using a 50 % threshold as the distinction between hail and hail-free events, the same as used in model construction.Validation was done for the HDT overall.The verification of the CM for cumulonimbus detection shows that of the 52 cumulonimbus events analyzed, the CM identified 48 correctly.Only four cases were considered nonconvective, none of which produced hail precipitation.In addition, only one cumulonimbus-free event (cirrostratus) out of 26 such events included in the verification database was erroneously classified as cumulonimbus.This was later correctly filtered by the HM.It can therefore be said that the CM does not filter out any hail event.
Table 8 shows average likelihoods of the two algorithms for the 26 events of each type included in the verification database.For the CM, greater likelihoods were obtained for the cumulonimbus events, with a likelihood near 100 % for hail events.The HM revealed great sensitivity, with a large difference between hail and hail-free cumulonimbus events.The two algorithms gave very low likelihoods for the cumulonimbus-free events.Once it was shown that none of the cumulonimbus-free events were classified as hail, the verification focused on convective episodes.
A contingency table was then built (Table 9) with results of the global application of the HDT for the 2011 cumulonimbus events.Note that only four events were wrongly classified as hail events.Of these, two corresponded to well-developed convective clouds generating intense rainfall (greater than 30 dBz radar reflectivity).Moreover, six hail events were wrongly classified as nonhail events.
Skill scores were computed using data from the contingency table to investigate different aspects of model validity.Table 10 shows a POD of 76.9 % and FAR of 16.7 %, both satisfactory values.Heidke's skill score (HSS) and true skill score (TSS) were also computed.For both indices, a value of 1 is considered a perfect forecast.
These results are slightly worse than those achieved by the radar-based hail-detection algorithms.However, the advantages of using satellite data instead of radar data make this approach valuable for monitoring hailstorms.

Application of the HDT: case study on 12 August 2011
The model was applied to the case study at 14:00 UTC on 12 August 2011.The synoptic situation showed a low pressure center at 500 hPa with strongly convergent flow at the surface.This center moved from the southwest of the Iberian Peninsula to the northeast.There was an isolated cold air mass in the upper tropospheric layers (−10 • C at 500 hPa), with strong warm and moist advection in lower layers.These meteorological conditions favored development of intense storms over the peninsula.In fact, the initial storms to the west-southwest of the peninsula were detected around noon.Shortly afterward, more storms developed over the Iberian and Central Mountain systems.Figure 5 shows areas with high likelihoods of accumulated hail as derived by radar during the day in the MEV.The Iberian range was the area most widely affected by hail.At 14:00 UTC there were several convective cells over the central peninsula, in different stages of development.In the north, nearly transparent cirrus and nimbostratus clouds were detected over the eastern coast.In the southwest of the peninsula there were stratus clouds (Fig. 2).The radar in the MEV showed a number of storm cells with high hail likelihood over the province of Teruel (Fig. 5).To the west of the Iberian Mountain system, lower likelihoods of hail were evident, but this might be a result of poor radar coverage because of the long distance .
The result of applying the CM and HM to the case study is shown in Fig. 6.The HDT results are shown on a probability scale, and are obtained by applying the HM to pixels identified as cumulonimbus by the CM (using a threshold of 50 %).Convective cells were detected with high hail likelihood in the southern Iberian Mountain system and around the Central Mountain system.These areas also agree with the radar for hailstorm presence (Fig. 7).

Conclusions
A daytime HDT was introduced for the summer months, using MSG data and applying two logistic regression models.The stepwise method was used to input variables for algorithm definition, with first-order interaction between the predictive variables.The CM identifies convective cloud pixels, whereas the HM discriminates pixels with hail precipitation from other convective cloud pixels.
The following conclusions can be drawn.
-The CM includes nine explanatory variables.Meteorological interpretation of these variables reveals that the "ingredients" required to discriminate cumulonimbus clouds from other cloud types are as follows: WV in the troposphere (assessed by the WV channels), thermodynamic phase of particles at cloud top and cloudtop temperature (assessed by the NIR channels and channel 8.7 µm), and cloud thickness and size of particles at cloud top (assessed by channel 0.8 µm and NIR channels).
-The HM includes four explanatory variables.Meteorological interpretation of these variables reveals that the "ingredients" required to distinguish hail areas within a cumulonimbus cloud are as follows: the presence of elevated OT (as shown by the visible channels), WV in the upper troposphere and updraft speed, based on small R e (derived from the NIR).
-Preliminary application of the CM is crucial to filter cloud tops formed by liquid water.The reason is that the HM is very sensitive to values from the NIR channels.Cloud tops formed by liquid water have small R e values and high reflectances in the NIR, so they would eventually be classified as hail.
-Analysis of the cloud physical properties indicated that cumulonimbus clouds are characterized by high LWP and OT values, plus varying R e values, depending on the stage of cell development.The occurrence of hail within cumulonimbus clouds is linked to lower R e values (30-50 µm) of particles at cloud top.Also, OT and LWP were not particularly useful parameters for identifying hail within a convective structure.
-Validation of the HDT, using independent data from 2011, gave a POD of 76.9 % and FAR of 16.7 %.
Although these results are slightly worse than those characterizing hail detection from ground-based radars, the HDT is recommended for application to areas with spatially limited radar coverage .The HDT will be used for tracking and nowcasting of hailstorms in real time, using radar data to implement the system.Finally, the effective spatial and temporal coverage of this tool will allow recording of hailstorms in the Iberian Peninsula in a database, for in-depth investigation of regional hail climatology.

Fig. 4 .Fig. 4 .
Fig. 4. Left: LWP (top), OT (middle) and Re(bottom) as in Fig. 3 for hail (black) and no hail (gray) pixels.Right: LWP (top), OT (middle) and Re(bottom), black lines correspond to areas with high likelihood of hail according to radar.Circled area shows radar range.26 Fig. 4. Left: LWP (top), OT (middle) and R e (bottom) as in Fig. 3 for hail (black) and no hail (gray) pixels.Right: LWP (top), OT (middle) and R e (bottom), black lines correspond to areas with high likelihood of hail according to radar.Circled area shows radar range.

Fig. 5 .
Fig. 5. 12 August 2011.Left: radar-based accumulated hail likelihood for entire day in MEV.Right: radar-based hail likelihood at 14:00 UTC.Convective cells with high hail likelihood are marked in red.Circled area shows the 100 km radar range

Table 2 .
Parameters selected by logistic regression for CM: metric explanatory variables input to the model, β coefficients used to extract the sign, and Wald parameter used to compute the weight.The symbol "•" represents multiplication between variables.

Table 3 .
Parameters of the global fit of the model for the CM.Chisquared test for −2LL reduction and several R 2 measures.

Table 4 .
Contingency table of database for CM.

Table 5 .
Parameters selected by logistic regression for HM: metric explanatory variables input to the model, β coefficients used to extract the sign, and Wald parameter used to compute the weight.The symbol "•" represents multiplication between variables.

Table 6 .
Parameters of global fit of the model for HM.Chi-squared test for −2LL reduction and several R 2 measures.

Table 7 .
Contingency table for HM database.

Table 8 .
Average results of algorithms for the 26 events of each type included in the verification database.CM shows cumulonimbus probability and HM hail precipitation probability.

Table 9 .
Contingency table for verification of HDT by cumulonimbus events in summer 2011.