Probabilistic forecasts of winter thunderstorms around Amsterdam Airport Schiphol

Abstract. The development and verification of a probabilistic forecast system for winter thunderstorms around Amsterdam Airport Schiphol is described. We have used Model Output Statistics (MOS) to develop the probabilistic forecast equations. The MOS system consists of 32 logistic regression equations, i.e. for two forecast periods (0–6 h and 6–12 h), four 90×80 km2 regions around Amsterdam Airport Schiphol, and four 6-h time periods. For the predictand quality-controlled Surveillance et Alerte Foudre par Interferometrie Radioelectrique (SAFIR) total lightning data were used. The potential predictors were calculated from postprocessed output of two numerical weather prediction (NWP) models – i.e. the High-Resolution Limited-Area Model (HIRLAM) and the European Centre for Medium-Range Weather Forecasts (ECMWF) model – and from an ensemble of advected lightning and radar data (0–6 h projections only). The predictors that are selected most often are the HIRLAM Boyden index, the square root of the ECMWF 3-h and 6-h convective precipitation sum, the HIRLAM convective available potential energy (CAPE) and two radar advection predictors. An objective verification was done, from which it can be concluded that the MOS system is skilful. The forecast system runs at the Royal Netherlands Meteorological Institute (KNMI) on an experimental basis, with the primary objective to warn aircraft pilots for potential aircraft induced lightning (AIL) risk during winter.


Introduction
In the Netherlands winter thunderstorms mostly form when cold polar air flows over the relatively warm North Sea, resulting in an unstable stratification.As a result deep convection may occur, potentially leading to thunderstorms.Little research has been carried out on winter thunderstorms in the Netherlands; a larger amount of research has been done in Japan.Kitagawa and Michimoto (1994) for instance describe the meteorological conditions and characteristics of Japanese winter thunderstorms.The most important feature is that thunderstorm clouds in winter are less deep than in summer, so the convective activity takes place at lower heights.
Winter thunderstorms are quite rare in the Netherlands, but still they can cause problems to aviation.Pilots have difficulties in estimating the risk of flying trough cumulonimbus clouds in winter, because these clouds are less deep than in summer and may seem harmless, but they could unex-Correspondence to: M. J. Schmeits (schmeits@knmi.nl)pectedly produce lightning.It is thought that about 90% of the discharges to aircrafts are initiated by the aircraft itself (Uman and Rakov, 2003).This is called aircraft induced lightning (AIL).When planes are struck by lightning, this may damage the plane which results in costly reparations or delays.That is why airports would like to have a warning system for the occurrence of lightning.There are several thunderstorm indices, like the Boyden index or CAPE (Table 1), which are quite good winter thunderstorm predictors when regarded individually.However, in this study thunderstorm indices are combined to create a probabilistic forecast system for winter thunderstorms.
In the United States, automated probabilistic forecasts for (severe) thunderstorms from 6 to 72 h in advance have been provided by the National Weather Service for a long period of time (e.g., Hughes, 2001).To our knowledge, in Europe a probabilistic forecast system like this only exists in the Netherlands.It forecasts (severe) thunderstorm probabilities in summer from 0 to 48 h in advance for 12 regions of 90×80 km 2 each (Schmeits et al., 2008).In this study a probabilistic 0-12 h forecast system for winter thunderstorms has been developed for four 90×80 km 2 regions Published by Copernicus Publications.dz square root of mean 6-h radar precipitation sum above 3 and 10 mm/h 7 around Amsterdam Airport Schiphol (Fig. 1b), including its 3 holding beacons as well (Slangen, 2008, hereafter referred to as S2008).
To develop the forecast system, statistical relationships have been determined between the occurrence of lightning and a set of potential predictors that were calculated from postprocessed output of the NWP models HIRLAM and ECMWF, and an ensemble of advected lightning and radar data (0-6 h projections only; Schmeits et al., 2008).
To determine the occurrence of lightning, a qualitycontrolled total lightning data set from the SAFIR network has been used.This dataset only contains those lightning discharges that match with radar echoes (S2008).An event is called a thunderstorm event if ≥1 lightning discharge is detected in a region in a 6 h period.
Model output statistics and logistic regression (Wilks, 2006) have been used to develop the forecast system.The run frequency of the MOS system is 4 times per day, and the forecast periods are 0-6 h (referred to as +6 h) and 6-12 h (referred to as +12 h).The equations are valid for the cold half-year only, i.e. from 16 October to 15 April.
In Sect. 2 the methodology is treated and in Sect. 3 we present the results.Finally, Sect. 4 contains some concluding remarks.

Logistic regression
The derivation of the MOS equations has been done using the method of logistic regression (Wilks, 2006).According to this method the probability Pr that an event y occurs is: . (1) The predictors x i (i=1, 2, ..., n) are selected via a so-called forward stepwise selection method and the regression coefficients a i are determined using the maximum likelihood method (Wilks, 2006).For an illustration of the logistic regression method the interested reader is referred to Fig. 8 of S2008.
Separate thunderstorm forecast equations have been derived for each region, projection and run time, resulting in a total of 32 forecast equations.No pooling has been performed in the derivation of the equations, except for the initial reduction in the number of potential predictors to assure some spatial consistency.Only data from the cold half-years from 16 October 2004 to 15 April 2007 have been used.The climatological probabilities of that period (see Fig. 11 of S2008) have been used as the reference forecasts in the verification of the regression equations (Sect.3.3).The climatological thunderstorm probabilities range from 2-4% in the 21:00-03:00 UTC period to 4-7% in the 15:00-21:00 UTC period, dependent on the region.

Potential predictors
As in Schmeits et al. (2008) several thunderstorm indicators from HIRLAM and from the ECMWF model, and an ensemble of advected radar and lightning data were used as potential predictors.The set not only included the traditional indices, like Jefferson or Boyden, but also other possible indicators of thunderstorms, like wind shear or temperature advection.
The HIRLAM predictor set consists of several traditional thunderstorm indices, like the Boyden index, the Jefferson index, CAPE and the Total Totals index.For each of the indices the minimum, average and maximum value was calculated using all grid points in each of the regions, and these values are used as potential predictors.This results in 51 potential predictors.The horizontal and vertical resolution of HIRLAM was changed in October 2006 (S2008).In each equation the latest HIRLAM cycle (i.e.00:00, 06:00, 12:00 or 18:00 UTC) is used.From the ECMWF model only the 12:00 UTC cycle is used.Various predictors, like vorticity or the square root of the 3-and 6-h convective precipitation sum, are included in this predictor set, which in total consists of 25 potential predictors.The horizontal and vertical resolution of the ECMWF model was changed in February 2006 (S2008).
Additionally to the NWP model predictors, observational predictors were added to the +6 h forecast potential predictor set.As in Schmeits et al. (2008), an 18 member ensemble of advected radar and lightning data was used.The lightning and radar images of 02:40, 08:40, 14:40 or 20:40 UTC were taken and advected 6 h ahead, using HIRLAM 700 hPa wind vectors and vectors calculated from previous radar images.The HIRLAM and radar vectors were each varied in length (25% shorter and 25% longer) and direction (+10 • and −10 • ), which resulted in a total ensemble size of 18 members.It was found in Schmeits et al. (2008) that most of the selected predictors were based on statistics of the 9-or 18member ensemble instead of individual ensemble members.Therefore, these statistics were used to create potential predictors.In this study a few new potential predictors were added, such as the square root of the 6-h radar precipitation sum above certain intensity thresholds (3, 10 and 30 mm/h).These predictors resemble the ECMWF convective precipitation predictors.15 Predictors of the radar and lightning advection ensemble were used as potential predictors.Besides the new predictors, the radar advection algorithm has been improved by the addition of a clutter filter that discriminates between rain and false radar echoes (S2008).

Selected predictors
Using forward stepwise selection (Wilks, 2006), 13 different predictors have been selected in the 32 forecast equations.However, many of these predictors show similarities.The most selected predictors are shown in Table 1.

Example case: 9 November 2007
To demonstrate the forecast system an example case is described.In the morning of 9 November 2007 small but heavy showers approached the Netherlands from the northwest.In Fig. 1a the 0-6 h probability forecast for 03:00-09:00 UTC is given.The probability is highest in the region WMS (92%).The other regions show less high probabilities, but still values above the climatological values (Sect.2.2).The observed number of lightning discharges is given in Fig. 1b.Lightning discharges were detected in WMS and MMS, while no discharges were detected in WMN and MMN.

Verification results
In this subsection objective verification results are shown for both forecast periods using the Brier skill score (BSS; Wilks, 2006).The BSS is computed for 30 different samples as follows.First, the data set is randomly divided into a dependent set ( 23 part of the data) and an independent set ( 1 3 part of the data), and this is done 30 times.As a lower limit for the number of cases with lightning occurrence in each independent set the amount of 4 cases is taken, to assure a certain level of statistical significance.The regression coefficients are calculated for the 30 dependent sets using the previously selected predictors.The coefficients differ slightly, because the sets are different.For every region and every time period, the BSS was calculated for the 30 independent sets.The Brier skill scores are shown in the boxplots of Fig. 2 and indicate whether the selected predictors can be used to forecast the occurrence of lightning.In Fig. 2a the boxplots are shown for the +6 h forecasts and in Fig. 2b for the +12 h forecasts, both for region MMN.Verification results for the other regions are similar and can be found in S2008.Thirty samples turned out to be sufficient, because only a small sampling error results when comparing two boxplots of each 30 samples corresponding to the same equation.
All the logistic regression equations generally have Brier skill scores above the 0-skill line (Fig. 2; Figs.24-27 of S2008).This clearly indicates that the forecast system developed here is skilful compared to climatology.The highest skills were found for the day and the lowest skills for the night.Another way to examine how skilful the equations are, is by testing them on a completely independent data set.This has been done for the data set from 16 October 2007 to 15 April 2008.From reliability diagrams based on these independent data, it can also be concluded that the forecast equations generally are skilful (Figs.29-32 of S2008).

Concluding remarks and outlook
Both in the winter (Table 1) and in the summer thunderstorm forecast system (Schmeits et al., 2008), the CAPE, the Boyden index and the convective precipitation predictors have often been selected.Only in the summer system, the Jefferson index and lightning advection predictors were often chosen.The Jefferson index is probably not selected in winter because it uses the temperature at 500 hPa, and winter thunderstorms generally do not reach that level.The lightning advection predictors are not selected in winter, probably because winter thunderstorms generally have shorter life times and less lightning activity than summer thunderstorms.This decreases the quality of a predictor that advects lightning activity, since it would probably overestimate lightning occurrence.
To determine the skill of the forecast system, we computed the Brier skill score for independent data sets.The Brier skill scores indicate that the equations generally have skill compared to climatology.The best skills are found for the day, less good but still generally positive skill scores for the night.
Possible future research could be to extend the MOS system to forecast projections of +24 h or even +48 h.Besides, the quality of the forecast system might be improved by taking the predictors over an area that is larger than the region for which the forecast is done.Finally, new potential predictors might be investigated, for instance from other observational sources like Meteosat Second Generation (MSG) and, in the future, from a short-term ensemble prediction system and from the non-hydrostatic NWP model HARMONIE.

Figure 1 .
Figure 1.(a) 0-6 h probability forecast (%) of ≥1 lightning discharge for 9 November 2007, 03:00-09:00 UTC.(b) Number of discharges per region in the same period.The location of Amsterdam Airport Schiphol is indicated by a cross in region MMN.The following letters are used for the regions(Schmeits et al., 2008): W stands for west, M for middle, E for east, X for extreme, N for north, and S for south.

Figure 2 .
Figure 2. Tukey boxplots showing the Brier skill scores of 30 independent samples as a function of the central verification time for region MMN, and for (a) +6 h and (b)+12 h forecast projection.The lower limit of the box gives the 0.25 quartile (q 0.25 ) and the upper limit gives the 0.75 quartile (q 0.75 ) of the data.The line inside the box is the median (q 0.50 ).The whiskers indicate the spread of the data, apart from outliers and extremes.Outliers are represented by an open circle (o), when they are 1.5-3 boxlengths away from q 0.25 or q 0.75 .A value is labelled extreme and marked with an asterisk (*) when it is more than 3 boxlengths away from the end of the box.

Table 1 .
Overview of the most selected predictors that are included in the indicated number of the 32 thunderstorm forecast equations.Here z is the (geopotential) height, T is the temperature ( • C), g is the acceleration due to gravity, LFC is the level of free convection, LNB is the level of neutral buoyancy, T v is the virtual temperature, and env. is short for environment.