Light (LGT) to moderate (MOD) aircraft icing (AI) is frequently reported at Cold Lake, Alberta, but forecasting AI has been a big challenge. The purpose of this study is to investigate and understand the weather conditions associated with AI based on observations in order to improve the icing forecast. To achieve this goal, Environment and Climate Change Canada in cooperation with the Department of National Defence deployed a number of ground-based instruments that include a microwave radiometer, a ceilometer, disdrometers, and conventional present weather sensors at the Cold Lake airport (CYOD). A number of pilot reports (PIREPs) of icing at Cold Lake during the 2016/17 winter period and associated observation data are examined. Most of the AI events were LGT (76%) followed by MOD (20%) and occurred during landing and takeoff at relatively warm temperatures. Two AI intensity algorithms have been tested based on an ice accumulation rate (IAR) assuming a cylindrical shape moving with airspeed υa of 60 and 89.4 m s−1, and the Canadian numerical weather prediction model forecasts. It was found that the algorithms IAR2 with υa = 89.4 m s−1 and IAR1 with υa = 60 m s−1 underestimated (overestimated) the LGT (MOD) icing events, respectively. The algorithm IAR2 with υa = 60 m s−1 appeared to be more suitable for forecasting LGT icing. Over all, the hit rate score was 0.33 for the 1200 UTC model run and 0.6 for 0000 UTC run for both algorithms, but based on the individual icing intensity scores, the IAR2 did better than IAR1 for forecasting LGT icing events.
Inflight icing occurs when an airfoil impacts with supercooled liquid drops that freeze on the surface, which can degrade the performance of an aircraft to the point of lost control. The type of icing, such as clear, rime, or mixed, depends on the cloud particle size distribution N(D) of the supercooled drops and temperature. The severity or intensity of icing is usually characterized as trace, light, moderate, and heavy, but the definition varies. The icing intensity definitions for pilot reports (PIREPs), which appear in the Aeronautical Information Manual (AIM), are very complex. They depend on many factors including the size and shape of the aircraft, ice accumulation rate, and duration of exposure to the icing environment. The Federal Aviation Administration (FAA) certification rule, based on the Code of Federal Regulations Title 14, Chapter 1, Part 25, Appendix C (FAR 25-C), only includes the smaller drops with a mean effective diameter less than 40 μm, and was developed for accretion of ice on a nonrotating cylinder moving at 89.4 m s−1 (174 kt or 200 mi h−1) with a diameter of 3 in. (Jeck 2001, and references therein). Now there is recognition of the importance of larger drops (>50 μm) usually referred to as supercooled large drops (SLD) as described in Cober et al. (2001, 2009) and Cober and Isaac (2012). Currently these finding are incorporated into the Appendix O icing envelope for SLD (https://www.gpo.gov/fdsys/pkg/FR-2014-11-04/pdf/2014-25789.pdf).
The icing intensity scale (IIS) that was originally adapted in the 1940s by the U.S. Weather Bureau [now incorporated into the National Oceanic and Atmospheric Administration (NOAA)] was defined as 0–1 g cm−2 h−2 [trace (TR)], 1–6 g cm−2 h−2 [light (LGT)], 6–12 g cm−2 h−2 [moderate (MOD)], or >12 g cm−2 h−2 [heavy or severe (SEV)]. The FAA Aviation Weather Research Program Inflight Icing Product Development Team proposed a new approach for determining icing intensity that depends on an alternative ice accumulation rate in centimeters per hour using several airfoil shapes including a cylinder (Politovich 2009). Using this method, the IIS is defined as 0.6–2.5 cm h−1 (light), 2.5–7.5 cm h−1 (moderate), and >7.5 cm h−1 (severe). This is the definition used in the FAA Pilot Guide: Flight in Icing Conditions (AC 91-74B, dated 8 October 2015) and hence, this IIS is adapted in this paper.
The implementation of these methods requires accurate prediction or measurements of cloud microphysical quantities, such as the median volume diameter (MVD), liquid water content (LWC), and other standard meteorological quantities. Although numerical weather prediction (NWP) models can be sophisticated enough to be able to predict these quantities, the accuracy of their prediction is still a big challenge, and hence methodologies were developed in the United States that incorporate more information such as satellite, surface station, and radar data for predicting icing potential (e.g., Bernstein et al. 2005). There have been improvements in the NWP models, but they are not adequately tested for their ability to predict icing intensity. Thompson et al. (2017), using the U.S. Weather Bureau ice intensity scale mentioned earlier and WRF Model data, found that their method overestimated the TR and SEV icing events and underestimated the LGT and MOD icing events when compared with PIREPs. For this study, pilot reports confirming icing events and freezing precipitation that occurred at Cold Lake, Alberta, Canada, during the 2016/17 winter periods are examined with data collected during the same period using a multifrequency microwave radiometer (MWR), a ceilometer, and other specialized instruments such as conventional weather sensors and disdrometers. The icing events are also examined using aircraft icing intensity estimated using the two ice accumulation models mentioned above with the relevant forecasts obtained from the Canadian High Resolution Regional Deterministic Prediction System (HRDPS). The main objectives of this study include 1) to test the ability of the HRDPS model to predict bulk microphysical quantities and associated mean particle size distribution N(D) as observed at ground level, 2) to test the ability of the MWR to detect reported icing events, and 3) to test the ability of the HRDPS model to forecast reported icing events. The sections are organized as follows. In section 2 a description of the study area and the instrumentation are given. A brief description of the HRDPS model is given in section 3. Comparisons of the observed bulk microphysical quantities at the ground level against the HRDPS model predictions are provided in section 4. In section 5, the description of the ice accumulation algorithms is given along with some case studies, and finally the data analysis results are discussed in section 6.
2. Study location and instrumentation
The study area, the Cold Lake Regional Airport (CYOD), is located in northeastern Alberta, Canada (541 m MSL; 54°N, 110°W). The location and climatology of the area are described in Boudala et al. (2017). The geographical locations of CYOD and the surrounding areas are shown in Fig. 1. Figure 1b shows a collocated permanent meteorological hourly observation (METAR) site used by the Department of National Defence (DND). To the west and northwest of the airport, there are a number of lakes (Marie, Ethel, Crane, and Hilda) including large lakes to the northeast of the site (Cold and Primrose). These lakes are known to contribute to the formation of precipitation and other weather phenomena when the flow is from a northeasterly direction. The west and southwest sides of the airport are surrounded by the Beaver River valley with an east–west orientation, which is also known to contribute to the formation of various weather conditions at the airport.
The meteorological instruments used in this study are shown in Fig. 2 and include a PARSIVEL disdrometer that measures particles size and velocity in the size range (0.2–20 mm) (Loffler-Mang and Lahak 2001; Loffler-Mang and Joss 2000), a Jenoptic CM15K ceilometer, and the Vaisala FS11P and PWD22 present weather sensors that measure visibility, precipitation intensity Pr and type, etc. The complete descriptions of the instruments are given elsewhere (Boudala et al. 2017, 2014). Here brief discussions of the MWR, the Meteorological Particle Spectrometer (MPS), and the Fog Monitoring Device (FMD) model (FM-120) are given.
The Radiometric MWR model MP-3000A uses 35 frequency channels to obtain vertical profiles of temperature T, relative humidity (RH), LWC, and absolute humidity ρυ at about 3-min intervals and vertical resolutions of 50 m from the surface to 500 m, 100 m from 500 m to 2 km, and 250 m from 2 to 10 km. The Radiometric radiometer uses statistical methods, with retrievals being performed using a standard back-propagation algorithm for training and a standard feed-forward network for profile determination, as described by Solheim et al. (1998). The training datasets can be historical rawinsonde or simulated data using the neural network retrieval method (Ware et al. 2013). Thus, the retrieval method is highly site specific, and hence the accuracy of these retrievals is not well known, particularly for the vertical profiles of LWC (Turner et al. 2007; Ebell et al. 2010). The objective of using the MWR in this study is to test if this instrument captures icing events reported by pilots and how the measurements compare with model data. The MWR is set to retrieve brightness temperature at the zenith, an elevation angle of 20°, and azimuthal angles of 180°, 225°, 270°, and 315°.
The MPS is an optical array probe that uses 64 linearly arranged photodiode elements to measure the shadow images of particles falling through a collimated laser beam (Baumgardner et al. 2002; Bringi et al. 2018). The shadow image that forms on the array is stored only if at least one of the diodes is shadowed by more than 50%. Knollenberg (1970) first introduced this measurement technique for measuring cloud N(D), and instruments were normally used onboard research aircraft but were later adapted by Baumgardner et al. (2002) for ground-based measurements. The MPS is mounted on a turntable with a wind vane to keep the sensor array oriented perpendicular to the average wind vector. This version of the instrument measures hydrometer number, velocity, and N(D). The instrument can measure N(D) at a resolution of 25, 50, or 100 μm with maximum diameters reaching 6.2 mm depending on the resolution. The instrument used in this study has a 50 μm resolution with a maximum size of 3.1 mm, and the first bin of the distribution was ignored because of measurement uncertainty. For the actual data analysis, the effective sampling area Aeff depends on the effective array width (EAW) and depth of field (DOF) given as
where dr is the probe resolution, n is the number of diodes (64), bin is the number of shadowed bins excluding the two end points of the diodes ranging from 1 to 62, is the wavelength, and is the diameter of a given particle. The sample volume is calculated as
where is the particle fall velocity and is the sampling interval. Droplet Measurement Technologies (DMT)’s FMD model FM-120 is a cloud-particle spectrometer designed for use during ground-based or tower-based studies (Spiegel et al. 2012). It measures cloud particle size and concentration in the size range 2–50 μm by drawing air through a sampling tube by means of a special pump at rate close to 15 m s−1. The probe has a sampling area of 0.24 mm2. The principle of operation is based on collecting scattered laser light in a forward direction (3.5°–12°) and by converting the scattered light intensity into particle size by means of light scattering principles. Using the three probes (FMD, MPS, and PARSIVEL) the full range of particle spectra from 2 μm to 20 mm can be measured and the bulk microphysical quantities such as LWC and MVD can be derived. As will be discussed later, the data can be used for model validation applications.
3. Model description
The model data used for this study is based on a 2.5 km pan-Canadian version of the Global Environmental Multiscale-Limited Area Model (GEM-LAM), which is now referred to as the High Resolution Deterministic Prediction System (HRDPS). A full description of HRDPS is given in Milbrandt et al. (2016), including the recent upgrades. The model has 60 levels, and the vertical resolution of the model varies from ~110 m near the surface to ~4 km at 40-km height depending on the atmospheric pressure, temperature, and humidity. The model was implemented with 48 h integration, run four times per day (at 0000, 0600, 1200, and 1800 UTC), and the output data has 5-min time resolution. The Canadian Meteorological Center (CMC) regional analysis data are used to provide the lateral boundary and initial conditions. The unified cloudiness–turbulence scheme, MoisTKE, is used for simulating the boundary layer clouds (Bélair et al. 2005; Mailhot and Bélair 2002). The subgrid cloud fraction and amount are determined based on conservative thermodynamic variables related to the water vapor saturation deficit (Bélair et al. 2005). Major modifications to the older version of GEM-LAM include the implementation of the two-moment version of the multimoment bulk microphysics scheme (Milbrandt and Yau 2005a,b) to calculate gridscale precipitation and the spectra of hydrometeors. Within this microphysics scheme, the hydrometeor N(D) of each category is by a gamma distribution function given as
where , μ, and λ are the intercept, shape, and slope parameters, respectively. Microphysical growth rates are formulated in two moments of the N(D), that is, particle number concentration NT and mixing ratio or cloud water content. For the liquid phase case, it is assumed that μ = 2, and knowing the LWC and for cloud drops or precipitation size particles, can be calculated as
where is a complete gamma function. The median volume diameter can be defined as
It can be shown that the median volume diameter can be given as
In HRDPS, only the number concentration of CCN is considered to determine the number of cloud droplets that form and hence is set at 80 cm−3.
4. Some examples of modeled and observed bulk microphysical quantities at ground level
The important bulk microphysical quantities that are necessary to be predicted accurately by NWP models for this study are LWC and MVD. Since reliable independent in situ measurements of these quantities are not available at cloud level other than the MWR for LWC, the ground level measurements of LWC and MVD were used to compare with the model predicted values at the same level. One case on 21 June 2017 during precipitation is discussed. The HRDPS uses a gamma size distribution in a form given in Eq. (4), with a μ = 2, for both cloud (c) and rain drops (r). The model predicts both the LWC and number concentration for rain and cloud drops separately. The modeled MVD was derived using Eq. (5), and the MVD was derived from observation data using the following procedures: since the MVD is the size of a droplet, below which 50% of the total water volume resides, it is derived using the following approach as suggested by DMT (DMT 2009). The accumulated mass, LWCn, is defined as
where is the size bin that varies from 1 to that satisfies Eq. (6), and LWC is the total water content. The MVD is computed by interpolating linearly between the size bins in a form
Figure 3 shows data obtained using several instruments and includes the hourly (METAR) and modeled data for 21 June 2017. The background light (BL) varied in the range 2–1000 cd m−2 representing normal night and day conditions (Boudala et al. 2012), and the visibility was reduced to near 4 km due to heavy rain during the night (Figs. 3a,e). The temperature varied from 5° to 20°C and the observed RH ranged from 50% to 100%, and the modeled predicted values reasonably agreed with observations (Fig. 3b). The observed maximum (minimum) LWCs for precipitation size particles were 1.5 (1.6) × 10−4 g m−3 and 2.1 (3.3) × 10−5 g m−3 for the MPS and PARSIVEL, respectively. No significant LWC was measured by the FMD (Fig. 3c). The model captured most of the later part (0500–1100 UTC) of the precipitation event (see Fig. 3e), but the predicted maximum (mean) LWC was 0.4 (1.5) × 10−5 g m−3, considerably lower than the observed values, and showed less variability (Fig. 3c). The maximum (minimum) MVD values measured using the MPS, PARSIVEL, and FMD were 2545 (75), 3612 (312), and 49 (3.5) μm, respectively. The model predicted MVD followed the MPS observation reasonably well with maximum (minimum) values of 1233.4 (562) μm (Figs. 3d,f). The maximum (minimum) concentrations measured using the MPS, PARSIVEL, and FMD were 1.77 × 105 (9.28) m−3, 4570 (0.81) m−3, and 4.44 (0.15) cm−3, respectively. The modeled maximum (minimum) concentrations of precipitation-sized particles were 2831 (0.13) m−3, which were significantly lower than the values observed by MPS, but relatively closer to the values observed by the PARSIVEL. No droplet size particles were predicted by the model on this day. The maximum concentration was measured by the FMD probe as would be expected, but the persistence of the concentration near 2 × 105 m−3 the entire day suggests that the probe was seeing aerosol particles with MVDs near 3–4 μm (Figs. 3d,f).
The maximum (mean) observed Pr values by the MPS, PARSIVEL, FS11P, and PWD22 were 17 (2.3), 18 (3.2), 22 (3.2), and 21.7 (2.9) mm h−1, respectively, and the precipitation type reported by the PWD22 was mainly rain (R) (Fig. 3e). Some drizzle was detected by the PARSIVEL (not shown here). As mentioned earlier, the model captured the later part of the precipitation event. However, the predicted maximum intensity was 7.5 mm h−1, which is lower than the observed values, but the mean value was 2.8 mm h−1, which is very close to the observed values owing to the fact that the modeled Pr is less variable as compared to the observations (Fig. 3e).
The mean and 95% confidence interval (CI) of the observed and forecasted bulk microphysical parameters are calculated using central limit theorem (Pek et al. 2017), which is valid for a sample size > 30. The mean and CI of LWC for the entire period (g m−3) varied as 0.27 ± 0.02, 0.24 ± 0.03, and 0.23 ± 0.01 for the PARSIVEL, MPS, and HRDPS, respectively. The predicted mean value is close to the observation, but the observed values are much more variable. The mean and 95% confidence interval for MVD (μm) varied as 26 ± 1, 359.6 ± 22.2, 1211.4 ± 23, and 1025.3 ± 13 for the FMD, MPS, PARSIVEL, and HRDPS, respectively. Similarly, the predicted MVD for precipitation-sized particles is close to the values measured by the PARSIVEL, but the model missed the smaller-sized particles. As expected the largest concentration is measured by the FMD ≈ 0.25 cm−3 (2.5 × 105 m−3) followed by the MPS (3.2 × 103 m−3), and 350 m−3 by the PARSIVEL. The averaged model predicted concentration was 0.001 cm−3, which is about 3 times lower than the measured values by the MPS, but the model captures more small particles than the PARSIVEL, which is consistent with the distribution given in Fig. 4. As mentioned earlier, since the observed values of the bulk microphysical parameters vary much more than the predicted parameters, on average the predicted parameters reasonably agreed with the mean observations within the time scale considered here, particularly for LWC and MVD.
Figure 4 shows some comparisons of hourly averaged predicted and measured N(D) using the FMD, MPS, and PARSIVEL probes and the associated observed mass fraction (MF) distribution for data collected on 21 June 2017. The figure shows the full N(D) (3 μm < D < 20 mm) starting from near the onset of the precipitation (Fig. 4a) to the end of the precipitation period (Fig. 4i; see Fig. 3). Generally, the MPS measures smaller drops than the PARSIVEL, and the PARSIVEL appears to underestimate the concentration at the smaller end of the spectrum (D < 600 μm), but there is a good overlap between the 600 and 3000 μm size range, which is consistent with the finding reported by Bringi et al. (2018). The maximum mass fraction measured by the PARSIVEL probe for particles D < 500 μm (drizzle size) was 5% as compared to 20% for the MPS probes. Contrary to the assumed gamma type N(D) in the model, the observed N(D) using the MPS probe suggests a nearly exponential N(D) with an increase in concentrations for D < 500 μm, and this is consistent with a recent study reported by Thurai et al. (2017). At the smaller end of the size spectra, the FMD observation suggests two peaks, one near 5 μm and one near 25 μm, and this is better illustrated in Fig. 5, which shows hourly averaged mass distributions. The model missed the onset of the precipitation events (Figs. 4a,b) and generally underestimated both the larger and smaller ends of the spectrum.
5. Aircraft icing
a. Icing intensity algorithms
where is the collection efficiency, LWC is the liquid water content (kg m−3), is the airspeed (m s−1), and is the density of ice assumed to be 800 kg m−3. The collection efficiency is calculated using the dimensionless Stokes number K and Langmuir’s parameter ϕ is defined as (Finstad et al. 1988)
where is the median volume diameter, is the density of water, is the diameter of the cylinder, is the dynamic viscosity of air, is the density of air, and is the Reynolds number. Following Finstad et al. (1988), the collection efficiency is defined as
where the coefficients A, B, and C are given in the following relationships:
Using the predicted LWC and from the HRDPS, the IAR values were calculated using Eqs. (8a) and (8b) for two values (89.4 and 60 m s−1). Traditionally m s−1 (200 mi h−1) was used for such calculations (Thompson et al. 2017; Jeck 2001), but through conversation with pilots at the 4Wing Canadian Air Force Base, at Cold Lake, Alberta, and based on PIREPs, most icing events occur during landing and takeoff when the speeds are generally lower than 89.4 m s−1. Consequently, the lower speed has been added for comparison. If the dynamic heating effect is not considered, which is the case here, the IAR increases with airspeed (not easily integrated in the model). The airframe surface temperature depends on many atmospheric factors including the aircraft speed as well as the ambient temperature. Cooling due to ice sublimation or warming due to deposition processes can occur (Mazin et al. 2001). Based on a nonrotating cylindrical model, the aerodynamic heating effect due to adiabatic compression of the air if the surface is free of ice, can reach approximately 1° and 3°C at an airspeed of 60 and 90 m s−1, respectively, but when the surface covered with ice the effect is relatively smaller (Mazin et al. 2001). An aircraft descending from an upper level with colder temperatures can get “cold soaked” and thus can be susceptible to icing at atmospheric temperatures well above 0°C.
b. Case studies
As described earlier in this study two icing intensity algorithms (IIA) based on IAR on a nonrotating cylinder with a diameter of 3 in. moving through clouds at 89.4 m s−1 as mentioned earlier and also 60 m s−1 have been tested by using the HRDPS forecasts of relevant parameters such as MVD, LWC, and other standard meteorological quantities as input. The models were then compared with the MWR and PIREPs. By virtue of the equations used in these two cylindrical models, the IIS TR, LGT, MOD, and SEV is given based on IAR. For the IAR1 model the IIS is defined as 0.6–2.5 cm h−1 (LGT), 2.5–7.5 cm h−1 (MOD), and >7.5 cm h−1 (SEV). For the IAR2 model, the IIS is defined as 0–1 g cm−2 h−1 (TR), 1–6 g cm−2 h−1 (LGT), 6–12 g cm−2 h−1 (MOD), > 12 g cm−2 h−1 (SEV). All available PIREPs, a total of 45, collected during the 2016/17 period were used in this study. As is typical, none of the PIREPs reported “no icing.”
1) 24 October 2016
On 24 October, light icing was reported by a C130 Hercules transport aircraft at 1849 UTC at about 217 m AGL, which is very close to the ground. Figure 6 shows the time series of modeled and observed vertical profiles of meteorological and bulk microphysical quantities. Both the model and MWR indicated enhanced LWC ranges of 0.01–0.36 and 0.01–0.6 g m−3, respectively, at lower levels (<1 km) and enhanced upward motion (w > 0 m s−1) aloft (Figs. 6a–c). The icing severity scale of TR, LGT, MOD, and SEV, based on the LWC threshold used by the U.S. Air Force (Thompson 1956), shown in Figs. 6a and 6b corresponds to green, yellow, dark red, and red colors, respectively. Based on the LWC threshold, LGT to MOD icing would be forecast depending on the time and height for both the model and MWR data (Figs. 6a,b). The height of the freezing level (0° isotherm) is near 1.5 km based on both modeled and MWR data, but the model showed a narrow-supercooled layer below about 1-km height consistent with the height that the icing event was reported. The modeled number concentration Nd and MVD ranges were 1–20 cm−3 and 1–60 μm, respectively (Figs. 6f,i). Both the MWR and HRDPS showed enhanced ρυ (~4 g m−3) and RH near the surface (Figs. 6g,h,j,k). The backscattered energy measured using the CHM15 ceilometer was totally attenuated, and it is consistent with a possible cloudy region identified using the predicted cloud fraction (CF) represented by 0.01 as shown in Fig. 6l delineating the supercooled cloud or fog.
Figure 7 shows the time series of the observed and HRDPS modeled Pr and type (Fig. 7a), temperature (Fig. 7b), humidity (Fig. 7c), wind speed (Fig. 7d), wind direction (Fig. 7e), observed background light (Fig. 6f), water vapor path (Fig. 7g), liquid water path (Fig. 7h), and visibility and ceiling (Fig. 7i). There was light rain and drizzle (Fig. 7a) during the fog event when visibility was reduced to less than 1 km. The prevailing wind directions during the fog and icing events were from east-southeast directions and the wind speeds were not very strong. However, there is a significant difference in the humidity measurements. Although the MWR uses a Rotronic RH/T sensor (MWR_rot) for surface RH and T measurements similar to the standalone sensor, the MWR measures smaller RH and warmer temperature which can affect the MWR profiling. Generally, the DewCell RH/T sensor can be taken as a standard measurement (Boudala et al. 2012) since it is well calibrated and kept in a Stevenson Screen. The lower MWR RH can partly explain why the humidity profile observed by the MWR is smaller than the expected humidity values in a saturated environment and the temperatures are relatively warmer. During the fog and icing events, enhanced liquid water path (LWP) was observed (Fig. 7h) and the fog event started during the nighttime (Fig. 7f,i), lifting near dawn, indicating that the fog may have been formed due to radiative cooling. Both the observed and model-predicted water vapor path (WVP) remained near 1 cm until 1600 UTC before the dissipation of the fog, but both increased sharply after 1600 UTC and during the icing event near 1848 UTC, particularly the data from the MWR that reached 1.8 cm as compared to the model predicted value of 1.5 cm (Fig. 7g). Similarly, both modeled and observed LWP increased starting from 0600 UTC, but the model predicted values were smaller by about a factor 2 (Fig. 7h).
Figure 8 shows the time series of IAR calculated assuming an airspeed of 89.4 m s−1 using the IAR1 and IAR2 algorithms described earlier for cloud drops (Figs. 8a,d) and raindrops (Figs. 8b,e), respectively, and cloud LWC and rainwater content (RWC; Figs. 8c,f). The IIS corresponding to the given model and LWC thresholds are also given similar to Fig. 6. The estimated icing intensity using IAR1 ranges from trace to heavy, but mostly moderate (Fig. 8a). Moderate icing is rarely reported in Cold Lake. Based on this data analysis, IAR1 seem to overestimate the intensity as compared to IAR2 (Fig. 8d) for a given aircraft speed, and cloud N(D). In this example, most of the particles contributing to the icing were MVD < 60 μm. The contribution of rain-size particles is relatively small (Figs. 8b,e,f). Based on discussion with the Air Force pilots at the base and many PIREPs, icing normally occurs during descent and take off when the values are relatively low (60–70 m s−1) for most aircraft; hence we have also tested the two models for an aircraft speed of 60 m s−1, and the results are given in Fig. 9. In this case all heavy icing events are removed from IAR1 and less moderate icing events were also forecast (Fig. 9a). If the IAR2 was used, however, all the moderate forecasts would be removed, and mostly light icing would be forecast (Fig. 9c), which would be much more reasonable for the Cold Lake weather conditions. The estimated icing intensity using IAR1 ranges from trace to heavy, but mostly moderate (Fig. 9a). Based on this data analysis, IAR1 seems to overestimate the intensity as compared to IAR2 (Fig. 9c) for a given aircraft speed and N(D).
Figure 10 shows vertical profiles of temperature (Fig. 10a), RH (Fig. 10b), ρυ (Fig. 10c), LWC (Fig. 10d), modeled updraft velocity w (Fig. 10e), potential temperature ϕ and the gradient with height (Fig. 10f), icing intensities (Fig. 10g), MVD (Fig. 10h), and cloud particle concentration (Fig. 10i) for the 24 October icing event (1849 UTC). The radiometer data shown here represents the data obtained at three angles: one zenith at an elevation angle of 90° and two azimuth angles of 270° and 315° at an elevation angle of ≈20°. Generally, both the observed and modeled temperature decreased with height above 1 km with the modeled data showing a sharp inversion between 0.5 and 1 km and also a minor inversion based on MWR data (270° and 315°), but not based on the Zenith data (Fig. 10a). Based on both the model and MWR data, the temperature of the layer where icing was reported varies between −3° and 2°C, which is relatively warm, but it should be noted that although that the icing was reported at this level, the accumulation could be initiated at a higher level where the temperatures are much colder (Fig. 10a) and where the MWR indicated LWC was present (Fig. 10d). Based on the model data, the atmosphere was saturated with respect to water at the level where the supercooled water was formed, but the RH sharply decreased with height reaching ~20% near 1.5 km (Fig. 10b). According to the MWR data, the RH remained constant (85%) up to about 1 km and decreased to 25% near 4 km (Fig. 10b). The ρυ varied from 0.5 to 5 g m−3 within the layer generally decreasing with height and characterized by positive vertical motion (Fig. 10e), but the model data shows some variation similar to the model RH (Fig. 10b). The observed LWC depends on the elevation and azimuthal angle and varied from 0.01 to 0.4 g m−3 in this example. The LWC measured at the 270° azimuthal direction was the highest. The modeled LWC reached 0.6 g m−3 but was limited to a narrow region. The values increased with height indicating a stable atmosphere, which is also confirmed by a strong positive vertical gradient (Fig. 10f). The model predicted icing intensity reached a maximum near 15 cm h−1 and 12 g m−2 h−1 for IAR1 and IAR2, respectively, at an altitude of 500 m where the temperature was near −3°C based on the model data (Fig. 10g). Depending on the model, this would imply heavy to moderate icing for IAR1 and IAR2, respectively. The model maximum droplet MVD and concentration reached 60 μm and 20 cm−3, respectively. The number of rain-size particles was very small.
Figure 11 shows comparison of the difference or bias between modeled and observed parameters using the MWR data shown in Fig. 10. The bias in LWC varied from −0.4 to 0.4 g m−3 depending on the height level and the MWR observation angle (Fig. 11a). Generally, the model underestimated (as compared to the MWR) the LWC below 4 km reaching a maximum of 0.4 g m−3 near the surface except between a narrow region within 0.4 and 0.6 km where the model significantly overestimated the LWC, particularly as compared to the zenith observation (Fig. 11a). The mean absolute error (MAE) [root-mean-square error (RMSE)] values of the bias in LWC for AGL < 4 km were 0.14 (0.18), 0.13 (0.15), and 0.11 (0.14) g m−3 for 270°, 315°, and zenith, respectively. The modeled temperatures were colder at lower heights (<1 km AGL) and warmer within heights between 1 and 2 km for all observation angles. The calculated bias using the observation obtained at zenith suggested a warmer layer for the model above 1 km AGL. The calculated MAE (RMSE) values of δT with respect to 270°, 315°, and 90° were 1.65°C (1.87°C), 1.45°C (1.68°C), and 1.9°C (2.2°C), respectively (Fig. 11b). The bias associated with relative humidity δRH was relatively larger; it varied from −40% at lower heights (1–1.8 km AGL), showing a drier layer for the model, reaching to about 50% within the layer between 2- and 3-km heights where the model is more humid. The calculated MAE (RMSE) values for δRH were 19% (22.26%), 19.4% (22.5%), and 18.5% (21%) for 270°, 315°, and 90°, respectively. The large difference in humidity could be attributed to uncertainties in both the model physics and MWR retrieval algorithms. The biases in the ρυ density vary from −1.5 to 1 g m−3 showing identical behavior as the RH profiles with relatively smaller MAE (RMSE) values of 0.55 (0.61), 0.52 (0.61) , and 0.49 (0.65) g m−3 for 270°, 315°, and 90°, respectively. The MWR data obtained at the 90° elevation angle gave slightly lower MAE values for δRH and δρ.
2) 26 October 2016
On 26 October, two icing events were reported. For the first one, light rime icing was reported by a CC130J Hercules cargo aircraft at 1442 UTC during descent between an altitude of 10 and 13 kft MSL (2.5–3.3 km AGL) and, for the second, engine icing was reported by a CF18 Hornet fighter jet aircraft at 1922 UTC during approach as the aircraft was passing through overcast cloud near an altitude of 186 m AGL. Figure 12 shows a similar plot as Fig. 6, but for the 26 October case. Both the model and MWR indicated enhanced LWC exceeding 0.25 g m−3 at some levels (<2 km), which implied MOD icing if the model temperature was used, particularly during the night between 0600 and 1000 UTC (Fig. 12l, low ceiling) at which time fog was reported. There was precipitation between 0200 and 500 UTC (see Fig. 13a), but the MWR did not indicate significantly enhanced LWC during this period. During the night (before 1400 UTC), both the MWR and the model have a warm layer (<2 km), but the model cold layer became wider during the day (Figs. 12d,e). The height of the modeled freezing level (0°C isotherm) sharply decreased to about 500 m until about two hours before the sunrise and then remained constant during the icing periods. The MWR data showed a similar trend, but the freezing level is rather variable probably because of the noise in the data and coarse vertical resolution. Generally, the model temperatures were colder and characterized by higher humidity values during the icing periods. The model-predicted clouds extended to higher levels during the night, as indicated by a bold pink thick line that delineates a possible cloudy region normally represented by 0.01 cloud fraction (Fig. 12l). This curve is not entirely visible at the surface level since the cloud fraction near the surface was 20% (not shown in the figure) indicating a very low cloud ceiling between 0000 and 1800 UTC. The backscattered ceilometer data suggests, however, that the cloud base heights are a little higher than predicted by the model, and it is also confirmed by the METAR data (thick black line; Fig. 12l). The predicted cloud top during the icing events was near 1.5 km for both events, but during the earlier event (1442 UTC), the pilot reported clear sky above about 3.3-km height, and the icing occurred between 2.5 and 3.4 km where the model suggested it was clear. The model, however, predicted a supercooled layer at a much higher level earlier before 1400 UTC. Contrary to the modeled data, the MWR data, suggests a supercooled layer extended to a much higher level, but the observed LWCs were relatively low (<0.1 g m−3). The modeled Nd and MVD ranges were 1–30 cm−3 and 1–49 μm, respectively (Figs. 12f,i).
Figure 13 shows surface conditions similar to Fig. 7, but for the 26 October case. The timing of observed precipitation between 0200 and 0500 UTC is well captured by the model (liquid phase), but the intensity is much smaller than observed (Fig. 13a). The surface temperature was well above freezing and decreased during the night and increased after sunrise (~1400 UTC), but the model predicted temperature increased during the night then followed the observations during the day when the icing events were reported (Fig. 13b). On the other hand, the RH was well forecast during and before the icing events including the foggy period (Fig. 13c). The MWR-rot sensor appeared to measure higher temperature and lower humidity on this day (not shown). The wind speed and wind direction measured with the WXT520 agreed quite well with the hourly METAR data (Figs. 13d,e). The wind speeds during the foggy events at night were well below 2 m s−1 when the prevailing wind direction were from east and southeast but increased later when the prevailing wind directions shifted to mostly westerly during the daytime when the icing events were reported. The model WVP derived from the dewpoint temperature and pressure forecasts agreed reasonably well as compared to the MWR observation, and both indicated enhanced WVP during the foggy events (Fig. 13g). The LWP measured with the MWR showed a stronger peak during the precipitation events indicating the influence of rain in the observed brightness temperature, but generally the model-predicted LWP followed the observation reasonably well (Fig. 13h).
Figure 14 shows the ice accumulation rates calculated assuming = 89.4 m s−1 similar to Fig. 8. The estimated icing intensity using IAR1 ranges from trace to heavy, but mostly moderate (2.5–7.5 cm h−1). As mentioned earlier, moderate icing is rarely reported in Cold Lake. Based on this data analysis, IAR1 seems to overestimate the intensity as compared to the IAR2 similar to the 24 October case. In this example, most of the particles contributing to the icing had MVDs smaller than 40 μm. The contribution of rain size particles is relatively small (Figs. 14b,e,f). Similar plots for υa = 60 m s−1 are given in Fig. 15. This example also demonstrates, as shown for 24 October, that less intense icing would be forecast at a lower aircraft speed, and if IAR2 was used, all the moderate icing forecasts would be removed, and mostly light icing would be forecast (Fig. 14c), which would be much more reasonable for the Cold Lake conditions.
Figure 16 shows the vertical profiles similar to Fig. 10, but for the icing event that occurred at 1442 UTC 26 October. Generally, the temperature decreased with height with a minor inversion near 0.6 km for the MWR and near 1 km just above the liquid layer for the model forecast. The freezing level according to the MWR data is located at a much higher level near 1.4 km as compared to about 0.6 km for the model data. The temperature where icing was reported varies approximately between −12° and −5°C based on both model and MWR data. There is a very little similarity between the RH measured and forecast (Fig. 16b). According to the MWR data RH remained constant (85%) up about 2 km and then decreased to 60%, but the model predicted higher RH near the surface increasing with height up to 200 m and then remaining constant up to 1 km and then decreasing to 75% at 3 km. The model forecast did not show icing conditions where the icing was reported since most of the model-based supercooled cloud layer was concentrated in the lower level. The MWR measured LWC up to about 0.05 g m−3 within the layer that icing was reported (Fig. 16d). The model-predicted icing intensity reached a maximum near 5.5 g m−2 h−1 and 7 cm h−1 for IAR1 and IAR2, respectively, at an altitude of 600 m where the temperature was near −2.5°C. Depending on the model, this would imply light to moderate icing for IAR1 and IAR2, respectively. The maximum droplet MVD and concentration reached 48 μm and 16 cm−3, respectively, and no significant rain-size particles were forecasted.
Figure 17 is similar to Fig. 11 but shows a comparison of the difference between modeled and observed parameters using the MWR shown in Fig. 16. The vertical profile of the bias in LWC has similar behavior as the case on 24 October but showed less variability ranging from −0.14 to 0.24 g m−3 with MAE (RMSE) of 0.08 (0.09), 0.07 (0.08), and 0.06 (0.08) g m−3 for 270°, 315°, and zenith observations (Fig. 17a). The model is colder for heights (<2 km) and warmer above this height based on the 90° and 315° observation angles. The calculated MAE (RMSE) values were 1.3°C (1.6°C), 1.2°C (1.4°C), and 2.4°C (2.7°C), for 270°, 315°, and zenith observations, respectively (Fig. 17b). The biases of the moistures variables RH and ρυ are also less variable as compared the first case. The MAE (RMSE) values for RH were 11.7% (13.7%), 12.7% (14%), and 13% (15%) and for ρυ 0.18 (0.22), 0.24 (0.3), and 0.24 (0.3) g m−3 for 270°, 315°, and zenith observations, respectively. The model has a moister layer at lower heights (<1.4 km) and a drier layer between 1.4 and 3 km depending on the direction of the MWR scan.
The second icing case that occurred at 1922 UTC was reported as engine icing near the ground during approach, and no icing intensity was reported. The temperature at this level was near zero based on both the model and MWR data. The RH remained constant near 90% up to about 2 km based on the MWR, and the model-predicted RH increased from 90% at the surface to 100% at a 1-km height. The maximum modeled LWC reached 0.4 g m−3 at about the 1-km height level. The corresponding maximum icing intensities, based on = 89.4 m s−1, were 6.5 g cm−2 h−1 and 8.5 cm h−1 for IAR2 and IAR1, respectively. These intensities are stronger than the first case, which would imply moderate to heavy icing, respectively. More reasonable intensity values can be calculated using = 60 m s−1 such as moderate to light icing for IAR1 and IAR2, respectively (see Fig. 15).
c. All icing cases combined
Figure 18 shows the combined 5-min modeled data for LWC plotted against the MVD for cloud drops for all cases including during times that icing events were reported. The Newton (1978) potential ice accumulation rate (IAR2) envelopes for two aircraft speeds 89.4 m s−1 (Fig. 18a) and 60 m s−1 (Fig. 18b), and similar envelopes for IAR1 (Figs. 18c,d) are also shown in the plot. The plot shows that the icing intensity increases with increasing LWC and sizes. This can be compared to Fig. 13 in Cober and Isaac (2012) that shows a similar plot but is based on in situ observation. Newton’s curves given in Fig. 18a, are identical to the ones given in Cober and Isaac. Visual examination of their plot indicates that the 30-s (approximately 3 km) averaged LWC values were lower than 0.5 g m−3, which would relatively agree with the MWR measurements. SEV icing is a rare phenomenon, and only one point was within the SEV icing envelope. On the other hand, as compared to the Cober and Isaac’s plot, the modeled LWC was underestimated at smaller MVD < 20 μm and overestimated for larger particle sizes reaching a maximum of 0.7 g m−3; hence there is a significant potential for forecasting SEV icing events, and the modeled light icing intensity for smaller size range droplets would be also underestimated. This may be partly attributed to the fact that on average the model underestimates the concentration of small particles, but this assumption needs further investigation. Furthermore, as shown in the figure, using a lower aircraft speed shifts the icing accumulation potential upward requiring higher LWC for icing to occur for the same MVD, particularly the severe icing case. As shown in Fig. 18c, in the IAR1, SEV icing can be forecast at much lower LWC of 0.2–0.3 g m−3 as compared to the 0.4–0.5 g m−3 for IAR2 at the same speed and MVD range (Figs. 18a,b). The two methods will be investigated using PIREPs in the following section.
d. Model and PIREP statistical comparisons
Using the 45 icing reports mentioned earlier and model forecasts at a 5-min time resolution, statistical tests were conducted based on two model initialization times (0000 and 1200 UTC) for comparisons. The nearest 5-min model forecast to the PIREP was selected for the verification.
Figure 19 shows the frequency occurrence of the icing events obtained from the PIREPs compared with models IAR1 and IAR2 for values of 89.4 and 60 m s−1 for all the 43 PIREP cases. Based on the PIREPS, the reported icing intensities were 4%, 76%, and 20% for TR, LGT, and MOD, respectively, and no severe cases were reported. The validations of the algorithms based on the 0000 UTC model run is given in Fig. 19a and the comparison for the 1200 UTC run is given in Fig. 19b. Generally, the events were reported during the day after 1400 UTC, so the 1200 UTC run is the closest to the events. As can be seen in the figure, about 76% of the PIREPs were LGT, but the IAR2 model with = 89 m s−1 for the 0000 UTC run forecast (54%) underestimated the event by about 29% and the 1200 UTC run forecast (28%) underestimated LGT icing events by 63% (Figs. 19a,b), and this is relatively consistent with the finding by Thompson et al. (2017). Using the lower speed of 60 m s −1 improves the light icing forecast, particularly for the 1200 UTC run forecast capturing close to 87% of the LGT icing events, and very few or no severe icing cases were reported depending on the forecast runtime (Figs. 19a,b). The IAR2 method significantly underestimated the MOD icing events for the 0000 UTC run forecast for both speeds (Fig. 19a), but when it was used with the higher speed the 1200 UTC run forecast overestimated the MOD icing events (54%; Fig. 19b). On the other hand, using the IAR1 algorithm with both air speeds significantly underestimated (overestimated) the LGT (MOD) icing events for the 0000 UTC run forecast, and also forecast SEV icing events that were not reported by the pilots. Thus, this study suggests that the ice accumulation rate based on the IAR1 model, particularly with = 89.4 m s−1 is not suitable for forecasting icing intensity using the HRDPS model for Cold Lake, but this could be different for different NWP models and locations. Generally, the model shows some skill forecasting icing, but the quality of the forecast depends on the model run time and assumed icing intensity algorithm.
Using the 45 PIREPs, statistical skill scores were calculated using a yes/no contingency table following Stanski et al. (1989). Since there are no “no” events in the PIREPs, only the hit rate (HR) or the probability of detection (POD) can be calculated for each IIS. For the overall yes events, the HR score was 0.33 for the 1200 UTC run and 0.6 for the 0000 UTC run, so the 0000 UTC run forecast appears to be better than the 1200 UTC run forecast particularly at the 60 m s−1 aircraft speed, although most of the icing events occurred close to the 1200 UTC run time. For the 0000 UTC forecast, the IAR1 (υa = 60 m s−1) method has the highest HR score of 0.33 for MOD icing events followed by 0.21 for LGT icing. The IAR2 the (υa = 60 m s−1) method showed a better HR score of 0.3 for LGT icing events. For the 1200 UTC forecast, only the IAR2 (υa = 60 m s−1) method showed some skill forecasting LGT icing with an HR score of 0.24. The overall HR scores mentioned above are similar to the ones reported by Guan et al. (2002) based on aircraft measurements and GEM model forecasts with three different cloud schemes. The work of Isaac et al. (2005) also showed the difficulty of verifying numerical model forecasts of icing with observations.
6. Summary and conclusions
Light (LGT) to moderate (MOD) aircraft icing (AI) is frequently reported at Cold Lake, Alberta, but forecasting the phenomena has been a big challenge. To understand the weather conditions associated with aircraft icing and improve the icing forecast, Environment and Climate Change Canada (ECCC) in cooperation with the Department of National Defence (DND) deployed a number of specialized ground based remote sensing instruments that included a multichannel microwave radiometer (MWR), CHM15k ceilometer, and other instruments such as the DMT Fog Measuring Device (FM-120) and a Meteorological Particle Spectrometer (MPS), a PARSIVEL disdrometer, and Vaisala present weather sensors at the CYOD airport in Cold Lake, Alberta. In this study, a number of pilot reports (PIREPs) confirming icing events that occurred at Cold Lake during the 2016/17 winter period and associated observation data were examined. The icing events were also examined using aircraft icing intensity estimated using two ice accumulation models based on a cylindrical shape approximation of an airfoil and the Canadian HRDPS model predicted bulk microphysical parameters such as LWC, droplet concentration Nd and median volume diameter (MVD), and temperature. These predicted bulk quantities were also compared with measurements at ground level during fog and precipitation. Two ice accumulation models (IAR1 and IAR2) with air speeds of 60 and 89.4 m s−1 were computed using the HRDPS forecasts and compared to against PIREPs. The major findings in this study include the following:
Surface observation data suggests that drizzle to precipitation size N(D) is better described by an exponential distribution function as opposed to the gamma type N(D) used in the HRDPS model.
At cloud level, the maximum LWC measured using the MWR within the layers that icing events were reported ranged from 0.04 to 0.4 g m−3 and the predicted LWC ranged from 0 to 0.7 g m−3. The MWR always gave some value of LWC during icing events, but the model did miss some cases.
The temperature and humidity of the cloud layers where icing was reported ranged from 0° to −23°C and from 50% to 100%, respectively.
Vertical motion predicted by the HRDPS model within the layers ranged from 0.01 to 0.075 m s−1 indicating upward motion as would be expected, but based both on the MWR and model data, the icing events occurred mainly in a stable atmosphere.
None of the reported icing events were formed due to melting snow followed by supercooling suggesting that the drop formation mechanism was through a condensation and collision–coalescence process; thus, local effects such as lakes and aerosol–cloud interaction could play a significant role.
Based on the PIREPS, the percentage of reported icing intensity events were 4%, 76%, and 20% for TR, LGT, and MOD, respectively, and no severe cases were reported.
The statistical validation of the HRDPS icing forecast showed that generally the model has some skill forecasting the icing events with an overall HR score of 0.33–0.6 depending on the runtime, but the quality of the forecast depends on the model run time and assumed icing intensity algorithm (IAR1 or IAR2). Based on this analysis, the IAR2 model with a lower speed of 60 m s−1 appears to be more suitable for forecasting LGT icing with a HR score 0.26 with some potential underestimation of MOD icing. The IAR1 algorithm with both air speeds significantly underestimated (overestimated) the LGT (MOD) icing events respectively for the 0000 UTC forecast, and also forecast SEV icing events, particularly with υa = 89.4 m s−1, thus making it not suitable for forecasting icing intensity in the Cold Lake environment. One of the reasons why the mode1 overestimated MOD icing, even predicting SEV icing under some conditions, could be associated with an overestimation of LWC.
It should be noted that these conclusions are only valid for the HRDPS model. The outcome could be different for different NWP models.
More studies using larger datasets with statistical performance tests are recommended. However, there is a problem associated with the lack of PIREPs that report “no” icing, because this can restrict the statistical tests available. If remote sensing techniques such as the MWR can be further validated, then perhaps the model performance for aircraft icing prediction could be validated using such remote sensing data.
This work was partially funded by the Department of National Defence (DND) and the Canadian National Search and Rescue New Initiative Fund (SAR-NIF) under SAR project SN201532. Mike Harwood and Robert Reed of ECCC helped with the installation of the instruments at Cold Lake. The authors also would like to extend their thanks to Raymond Dooley, Amy Slade-Campbell, and Gordon Lee of DND for providing the PIREPs, and also Randy Blackwell of DND for helping during the instrument installation process.