Forecasting severe ice storms using numerical weather prediction : the March 2010 Newfoundland event

The northeast coast of North America is frequently hit by severe ice storms. These freezing rain events can produce large ice accretions that damage structures, frequently power transmission and distribution infrastructure. For this reason, it is highly desirable to model and forecast such icing events, so that the consequent damages can be prevented or mitigated. The case study presented in this paper focuses on the March 2010 ice storm event that took place in eastern Newfoundland. We apply a combination of a numerical weather prediction model and an ice accretion algorithm to simulate a forecast of this event. The main goals of this study are to compare the simulated meteorological variables to observations, and to assess the ability of the model to accurately predict the ice accretion load for different forecast horizons. The duration and timing of the freezing rain event that occurred between the night of 4 March and the morning of 6 March was simulated well in all model runs. The total precipitation amounts in the model, however, differed by up to a factor of two from the observations. The accuracy of the model air temperature strongly depended on the forecast horizon, but it was acceptable for all simulation runs. The simulated accretion loads were also compared to the design values for power delivery structures in the region. The results indicated that the simulated values exceeded design criteria in the areas of reported damage and power outages.


Introduction
The Canadian Maritime Provinces and the northeastern states of the USA are often affected by severe ice storms caused by freezing rain.In 1998, one of the most severe and Correspondence to: P. Musilek (petr.musilek@ualberta.ca)extensive ice storms devastated densely populated areas in eastern Ontario and southern Quebec, Canada, and rural areas in northern New York and New England, United States.In March 2010, another severe ice storm struck the island of Newfoundland, Canada, leaving many people without power for several days.During this event, the most intense freezing rain occurred overnight on 4 March and throughout 5 March.The storm caused extensive power outages affecting about 7000 customers in 32 communities on the Bonavista Peninsula and 22 in the northern Avalon Peninsula.The outages were caused by physical damage to about 250 power distribution lines and support structures that broke under the weight of accreted ice.Media reported that conductors were covered with ice up to 12 cm thick (Truro Daily News, 2010).
These two cases show the potential magnitude of freezing rain storms, and illustrate their technical, economic and societal consequences.It would be highly desirable to model and forecast such icing events, so that the consequent damages can be prevented or mitigated.The goal of this paper is to evaluate the forecasting capability of our Ice Accretion Forecasting System (IAFS), which is based on a Numerical Weather Prediction (NWP) model and an ice accretion model augmented with an intelligent freezing rain detection algorithm (Musilek et al., 2009).
Unlike many meteorological parameters, icing measurements are not standard at meteorological observation stations.Furthermore, there is still a lack of a reliable icing sensor that can work under all expected climatic conditions (Fikke, 2009).As a result, icing can rarely be objectively monitored.Similarly, forecasting models can seldom be quantitatively validated due to the lack of observational data.Nevertheless, a number of ice accretion models have been developed that can predict ice accumulations based on meteorological variables.Some of these ice accretion models, e.g. the morphogenetic or random walk scheme by Szilder (1993), are more suited for detailed icing studies rather than for operational use.Diagnostic models are usually based on the mass flux of water droplets reaching the surface of the icing object and the efficiency of the accretion process, e.g.Makkonen (2000).Modelling freezing rain accretion using diagnostic models involves the adoption of several simplifying assumptions that make their operational use less computationally demanding.First, freezing raindrops are assumed to be large enough so that their collision efficiency can be considered unity.Second, the wind is assumed to be perpendicular to the line, which is generally considered to produce the highest ice loads.Third, the accretion is assumed to be a circular cylinder.Thin ice accretions on overhead power line conductors are quasi-circular, as are large ice accretions on torsionally weak conductors.However, when the ice accretion is dry and the conductor is torsionally stiff, the ice forms on the upper, windward side, giving rise to a quasi-elliptical shape.At air temperatures within a few degrees of freezing, icicles tend to form.Because icicles intercept additional droplets and provide an opportunity for unfrozen surface liquid to freeze instead of being shed, ice loads with icicles can exceed those computed with a simple, cylindrical model.Fourth, melting is ignored.Finally, using all these assumptions, a simple model can be formulated to calculate the radial equivalent ice accretion thickness on a cylinder, which does not depend on the original size of the object (Jones, 1998).A survey of various approaches to ice accretion modelling can be found in Makkonen (1998).The author concludes that simple models are sufficient for freezing rain simulation, as long as the ice growth remains in the dry regime.
In order to forecast ice accretion due to freezing rain, future values of meteorological variables pertinent to glaze ice accretion must be known, including the freezing rain precipitation rate, wind velocity and air temperature.These values can be forecast (or simulated) with a NWP model.However, NWP models do not typically directly output the freezing rain precipitation rate.Consequently, we augment the NWP model with an intelligent freezing rain detection algorithm.This algorithm allows us to invoke the icing model only when the combined NWP model and intelligent algorithm indicate that freezing rain is actually falling.The idea of combining a NWP model with an icing model was suggested more than a decade ago (Vassbo et al., 1998).Due to the advanced microphysics schemes invoked in the model, the problem generally requires a fine spatial resolution, resulting in high computational costs.However, the rapid development of computer power during the last decade has made it possible to combine NWP and icing models, while the accuracy of NWP systems has significantly improved (Thompson et al., 2009).
NWP models were recently used to predict rime icing events in mountainous regions (Nygaard, 2009;Dierer et al., 2009).A NWP model, coupled with a universal mass-based icing scheme that resolves both rime and glaze icing was applied in northeast Bulgaria (Nygaard and Nikolov, 2009).
While the latter model captured the dominating weather situation that produced freezing rain well, it failed to localize the event due to anomalously high surface temperatures.The authors suggest that this error may have been caused by too strong vertical mixing in the planetary boundary layer (PBL) scheme, and poor representation of the surface stable layer in the model.
In this study, we apply our IAFS to simulate a forecast of the March 2010 Newfoundland ice storm, in order to evaluate its ability to capture and quantify such events.Although this simulation was carried out after the storm occurred, it was initialized and provided with boundary conditions from global model forecasts available up to 84 h prior to the main icing event.This paper is organized in six sections.Section 2 briefly describes the ice storm, while Sect. 3 provides an overview of the forecasting system.NWP simulations are compared to observations in Sect.4, and Sect. 5 provides a thorough analysis of the ice accretion estimates.Major conclusions and possible directions for future work are outlined in Sect.6.

Ice storm description
A surface low pressure system formed on 2 March 2010, over the Gulf of Mexico; it started to move towards the northeast, driven by the flow in the middle troposphere.The depression deepened on 3 March, while passing Cape Hatteras on the US east coast.The minimum sea-level pressure dropped from 1002 hPa on 2 March to 986 hPa on 3 March.The storm hit southeastern Newfoundland on 5 March, and brought rain and freezing rain to the exposed northeastern coast of the island.The slowly moving storm produced rain for about two days.The highest total precipitation exceeded 110 mm on the Avalon Peninsula.
In terms of the duration and maximum precipitation, this storm is comparable with the Great Ice Storm of January 1998.Another common characteristic was that both storms were slow moving.The major differences were the spatial extent of the icing, and the fact that, luckily, much of the area affected by the more recent storm has a relatively low population density.However, low population density can cause long restoration times due to the prioritization of line repair operations (Short, 2004).

NWP model setup
Weather simulations for the current study were performed using the WRF ARW model version 3.2.WRF is an open source, mesoscale, non-hydrostatic NWP model, developed for research and operational weather forecasting (Skamarock et al., 2008).Three nested model domains were used with grid sizes of 10.8 km, 3.6 km and 1.2 km, and grid dimensions of 47 × 52, 64 × 79, 118 × 163, respectively.The innermost domain, covering an area of 141.6 × 195.6 km, is shown in Fig. 1, along with the main power line on the Bonavista Peninsula.In order to reproduce the conditions of a real forecast, initial and boundary conditions were obtained from the North American Model (NAM) data products, which are based on global model forecasts.To determine how far ahead the ice storm could be predicted, forecast horizons of 64, 52, 40, 28 and 16 h relative to the main event were considered.No observation nudging or data assimilation was performed.
For modelling icing on structures in the surface layer of the atmosphere using an NWP model, the choice of physical parameterizations is crucial.Since the main inputs to ice accretion models are freezing rain precipitation rate, liquid water content and wind speed, the most important parameterizations are the microphysics (MP), planetary boundary layer (PBL) and surface layer (SL).The choice of other schemes has a negligible impact on the simulation results.The long-wave radiation scheme, Rapid Radiative Transfer Model (Mlaver et al., 1997), and the surface scheme, Pleim-Xiu Land Surface Model (Xiu and Pleim, 2001), were used in the simulations.Convection is assumed to be explicitly resolved by the finest model grid, and thus no cumulus scheme was applied for the innermost domain, while the outer domains used the Grell-Dévényi parameterization (Grell and Dévényi, 2002).
The cloud microphysics was parameterized using Thompson's scheme (Thompson et al., 2004).This scheme is used frequently in icing modelling (Nygaard, 2009;Dierer et al., 2009).It is the best-tested, double-moment microphysics scheme, and the only double-moment scheme included in the operational core of the WRF model (NMM).Since the amount of liquid water in the model atmosphere is very sensitive to the microphysics parameterization, a singlemoment scheme would not provide sufficient accuracy.In addition to the mass concentration, the double-moment scheme also provides the number concentration of liquid and/or solid precipitation particles.
The boundary and surface layers were parameterized using the Quasi-Normal Scale Elimination (QNSE) scheme (Sukoriansky and Galperin, 2008), which is an alternative to Reynolds stress turbulence models.The method employs a successive scale elimination and calculates corrections to the viscosity and diffusivity for a given scale.It accounts for anisotropy in the turbulence caused by thermal stratification.It also resolves the stable surface layer, which is often present during freezing rain events, better than classical methods (Mellor and Yamada, 1982;Hong et al., 2006).

Icing model
The freezing rain model applied in this study is based on the Simple Model (SM) of Jones (1998).It calculates the radial equivalent ice thickness, R eq in [mm], which is the thickness of a radially uniform ice accretion on a conductor, with the same mass as the actual ice accretion.Because of the large size of raindrops, the radial equivalent ice thickness of freezing rain is independent of the conductor diameter.The total radial equivalent ice thickness is computed for an entire icing event as follows  the Simple Model to the METAR data (regular reports from airport meteorological stations), triggering the icing algorithm with actual observations of precipitation type.However, when forecasting icing with a NWP model, observations cannot be used to distinguish precipitation type, as they are not available at the time of forecast.Instead, the forecasting system has to rely only on information provided by the model.NWP models offer enough information to assess the precipitation type, namely the fraction of frozen precipitation and air temperature.Fraction of frozen precipitation is a dimensionless variable defined as the ratio of frozen precipitation (in the form of ice, such as snow and ice pellets) and total precipitation.To distinguish between rain and freezing rain, air temperature and fraction of frozen precipitation can be used as described in Ramer (1993).However, as shown by Musilek et al. (2009), diagnosis of precipitation type with fixed thresholds for these variables may not lead to satisfactory results.Consequently, these authors developed an IAFS (Ice Accretion Forecasting System) that uses a fuzzy algorithm to control engagement of the ice accretion model, only at times when freezing rain is identified by the algorithm.This approach was refined by Pytlak et al. (2010), using ASOS observations and a genetic algorithm to identify the optimal setup for the engagement function, based on model surface temperature, and the fraction of frozen precipitation.This advanced IAFS was incorporated into the present study.

Comparison of NWP outputs with surface observations
Before analysing the freezing rain ice accretion model, the relevant simulated meteorological parameters (precipitation, wind speed and temperature) were compared with surface measurements inside the innermost domain.This domain encompasses five meteorological stations operated by Environment Canada.Most of the stations only make daily observations.Only one station, Bonavista, provides data with hourly resolution, suitable for temperature and wind speed verification.All model values used for the comparison were transposed using bilinear interpolation from the original gridded WRF output to the site of the Bonavista station.The amount of freezing rain is the most important factor determining the intensity of glaze icing.However, since observations of freezing rain were not available at any of the five stations, we instead compared daily totals of precipitation observed on 5 and 6 March, with the total amounts provided by the NWP model.Inasmuch as the total precipitation may include rain, freezing rain, ice pellets and snow, this is not an ideal comparison.However, it may indicate whether or not the model has a systematic tendency to underestimate or overestimate precipitation in general.The amounts of total precipitation observed on 5 and 6 March are compared to the values provided by the NWP model in Tables 1 and 2. The tables also evaluate the measured and simulated amounts in terms of Mean Absolute Error (MAE), across all meteorological stations.
At Environment Canada stations, the climate day at first order or primary observing sites is defined by the 24-h period ending at 06:00 UTC.However, at volunteer observing sites, the standard rain gauges are read once a day at approximately 8 a.m.local time, which is 11:30 UTC.For consistency, the simulated daily precipitation totals were taken either at 06:00 UTC or at 11:30 UTC depending on the type of station.The largest errors appear at the Terra Nova station, which reported the highest amount of precipitation on both days.The model's tendency to underestimate total precipitation could be caused by the complex terrain surrounding the station, which is not well resolved at the spatial resolution we used.The MAE generally decreases with a shorter forecasting horizon, with one insignificant exception, namely the last simulation for 5 March.The model precipitation field (Fig. 2) appears to be flattened compared to the observational field, with smaller observed values overestimated by the model and vice versa.Hence, it appears that the model precipitation is underestimated at the Nat.Hazards Earth Syst.Sci., 11, 587-595, 2011 www.nat-hazards-earth-syst-sci.net/11/587/2011/The meteorological station in Bonavista belongs to the network of synoptic stations, and thus precipitation data are available at 6-h intervals.The comparison of the model and observations is shown in Fig. 3.The data are presented as cumulative values starting at 00:00 UTC on 4 March.The precipitation graphs for model runs started later than this are initialized with precipitation values observed at their onset.The simulation started closest to the beginning of the icing event is significantly more accurate than the earlier ones.On average, the total model precipitation was about 60% of the observed total precipitation.
For wind speed and air temperature, the only meteorological station that provides hourly data is Bonavista.However, the anemometer was apparently frozen during the main icing event, and so the wind speed data are missing during that period, preventing verification.The modelled 10-m wind speed at the location of the Bonavista weather station increased from about 11 m s −1 to about 16 m s −1 at the height of the storm and then declined to near 15 m s −1 at the end www.nat-hazards-earth-syst-sci.net/11/587/2011/Nat.Hazards Earth Syst.Sci., 11, 587-595, 2011 of the freezing rain period.Peak wind speeds as high as 23 m s −1 in the model occurred offshore and along Trinity Bay, Placentia Bay and Conception Bay.The wind direction was predominantly NNW, coinciding with the longitudinal axes of these bays.Since the maximum model precipitation rate was only a few mm h −1 , the strong winds made a very substantial contribution to the overall icing rate in the model.They may also have contributed significantly to the overall load on the lines, although we did not attempt to calculate the wind load.Air temperature observations are available, and their comparison with modelled values from two simulation runs can be found in Fig. 4. The first run was initialized at 00:00 UTC on 4 March, well before the event started.The resulting simulated temperatures are significantly lower than observed in the second part of the event, and the peaks in the time series are also shifted by about 6-10 h.The simulation initialized at 00:00 UTC on 5 March performed much better, with no time shifts and errors of 0.3 • C or less.

Forecast of freezing rain
The Simple Model was applied throughout the domain to calculate the glaze ice accumulation due to freezing rain.The IAFS algorithm (Musilek et al., 2009) was employed to discount precipitation that fell in frozen form.The best performing setup, IAFS5GA (Pytlak et al., 2010), was used for the simulations.This algorithm uses fraction of frozen precipitation and surface wet bulb temperature to determine an engagement function for the Simple Model.Both values are derived from the outputs of the WRF model.The distribution of the fraction of frozen precipitation over the simulation domain is shown in Fig. 5.The engagement function determines the freezing rain portion of the ice load, calculated with the Simple Model and using the total model precipitation.The results of the simulation initiated at 00:00 UTC on 5 March are shown in Fig. 6.The most affected areas (the Bonavista Peninsula and Conception Bay North) are well captured by the model.The maximum model radial equivalent ice thickness exceeds 33 mm.This value is lower than values reported in the media (up to 60 mm), but it is important to note that the total precipitation is underestimated by a similar factor.It is also possible that values reported in the media are exaggerated because of icicle fomation and ice accretion asymmetry.
This affected area of eastern Newfoundland belongs to the severe icing category according to the overhead systems design standard (Canadian Standards Association, 2006), and the lines in this category are supposed to be designed for 19 mm of radial equivalent ice thickness.The spatial extent of the area with simulated radial equivalent ice thickness exceeding 19 mm is overlaid on a map of Newfoundland municipalities in Fig. 7.All three regions with reported heavy damage on power lines -the Bonavista Peninsula, the northern Avalon Peninsula and Conception Bay Northcontain simulated ice loads exceeding the design threshold.Figure 8 shows the radial equivalent ice thickness along the main transmission line on the Bonavista Peninsula (the position of the line is shown in Fig. 1).The values at the beginning of the profile, corresponding to the northern part of the line, do not change because the line follows the coast at relatively low elevations.Further along, the line crosses an exposed area, which corresponds to the peak of the simulated icing load.The icing load then decreases towards the interior of the island.The largest icing loads are simulated for the longest forecast horizon, initialized on 3 March.The run initialized at 00:00 UTC on 4 March (forecast horizon of approximately 40 h relative to the main event) produces smaller loads than the other runs.However, all forecasts would have been useful in alerting the power utilities about the approaching storm, because they predict maximum ice loads approaching or exceeding the design limit.
Time series of the radial equivalent ice thickness at the point with maximal values are shown in Fig. 9.The simulation initialized on 3 March at 00:00 UTC produced an earlier onset of the icing compared to the other simulations,   by about 10 h.The simulations with later starts agree better with the reports that freezing rain was observed throughout 5 March.
The timing and duration of the storm are consistent with reports from Environment Canada (www.ns.ec.gc.ca, 2010).For the Bonavista peninsula, the report states that, "The precipitation fell mainly as freezing rain overnight Thursday and all day Friday".This time period corresponds, approximately, to the time window delineated by two vertical lines in Fig. 9.

Conclusions
In this paper, an Ice Accretion Forecasting System (IAFS), consisting of a NWP model and an ice accretion algorithm, was used to evaluate our ability to forecast the severe ice storm that took place in Newfoundland on 4-6 March 2010.The forecasting system consisted of the WRF model version 3.2, the Simple Model for ice accretion, and a freezing rain engagement function optimized by a genetic algorithm.Although the IAFS was used to simulate the event after it occurred, it was employed in a way similar to how it would be used as a forecasting tool.All data used to initialize the simulations was available as global forecasts up to 64 h before the onset of the event.The boundary conditions for the outer domain were also supplied by the global forecast data.
Meteorological variables simulated by the WRF model were compared with available surface observations.Although the precipitation timing was generally accurate, the total precipitation at some stations was underestimated or overestimated by as much as a factor of two, possibly because of erroneous estimation of the locations of precipitation maxima.In general, the total precipitation at locations with higher observed accumulations was underestimated.The wind speed was not compared with observations due to missing observed data.The accuracy of the surface air temperature simulation was acceptable for all model runs; it improved significantly for shorter forecast horizons.Ice storm timing and duration were properly simulated, between the night of 4 March and the morning of 6 March, in all model runs started after 3 March 00:00 UTC.The radial equivalent ice accretion thickness appears to be underestimated or overestimated by a factor similar to that for total precipitation.
In summary, based on our simulation of this very severe icing event, the IAFS appears to be a promising tool for operational forecasting of severe freezing rain events.Nevertheless, there is still room for improvement of the freezing rain precipitation forecast, which has a crucial impact on the modelled icing load.

Fig. 1 .
Fig. 1.The inner domain of the WRF model (141.6×195.6 km) and the main power distribution line on the Bonavista peninsula (red).
where index j denotes values in the j -th hour, P is the freezing precipitation rate [mm h −1 ], ρ 0 is the density of water [kg m −3 ], ρ i is the density of ice [kg m −3 ], U is the component of wind velocity normal to the line [m s −1 ], W is the precipitation liquid water content [kg m −3 ], and N is the duration of the ice storm[hours].The precipitation liquid water content is computed as a function of the precipitation rate according toBest (1950) W = 6.7 × 10 −5 P 0.846 .(2)The2002 ice storm in the southern United States was analyzed byJones et al. (2004).The authors applied www.nat-hazards-earth-syst-sci.net/11/587/2011/Nat.Hazards Earth Syst.Sci., 11, 587-595, 2011

Fig. 4 .
Fig. 4. Observed surface air temperature at Bonavista and temperature modelled with WRF at 2 m a.g.l. for simulations started on 4 and 5 March 00:00 UTC.

Table 1 .
Daily total precipitation (mm)observed on 5 March 2010, and simulated total precipitation for the same period.

Table 2 .
Daily total precipitation (mm)observed on 6 March 2010, and simulated total precipitation for the same period.