## Abstract

Short-term automated forecasts (nowcasts) of liquid water equivalent (LWE) values can be used to assist aviation deicing decision-making activities. Such decisions can mitigate hazards that cause losses of life and property and increase costs because of travel delays. The Weather Support to Deicing Decision Making (WSDDM) system provides LWE nowcasts and is currently deployed at several major airports in the United States. WSDDM produces these nowcasts in two steps. First, an equation relating radar reflectivity to LWE rate is calibrated by correlating radar and surface observations. Then, nowcasts of reflectivity are converted to nowcasts of LWE using this calibrated equation. This paper shows that the incorporation of the Dynamic and Adaptive Radar Tracking of Storms (DARTS) radar–based nowcasting method into WSDDM can provide more accurate and efficient nowcasts of LWE relative to the correlation-based nowcasting method currently used. Results of an evaluation considering approximately 92 h of data collected during four winter weather events show the incorporation of DARTS into WSDDM provides an approximate 14% average improvement in the accuracy of 60-min LWE nowcasts and reduces runtime by two orders of magnitude.

## 1. Introduction

The operation of aircraft during winter precipitation conditions is a significant safety concern. A relatively thin coating of ice on an aircraft wing can cause substantial loss of lift and increase in drag, making safe flight impossible (Rasmussen et al. 2000). Thus, the formation of ice must be removed or prevented using a deicing and anti-icing fluid prior to takeoff. Additionally, runways are often treated with an anti-icing fluid to prevent the accumulation of snow that can hamper safe and efficient airfield operations during winter storms (Rasmussen et al. 2003).

A key quantity used to assess the type and application frequency of deicing and anti-icing fluids is the liquid water equivalent (LWE) value. Accurate short-term forecasts (nowcasts) of LWE values are useful to facilitate efficient resource allocation and planning of airport deicing operations. The Society of Automotive Engineers’ (SAE) Ground Deicing Committee publishes standards for the four different types of aviation deicing fluids (Society of Automotive Engineers 2011a, 2011b). Based on these standards, the Federal Aviation Administration (FAA) publishes annual revisions to the official holdover time tables for all approved deicing fluids (Federal Aviation Administration 2011). The holdover time tables indicate the amount of time a deicing fluid can prevent ice formation on an aircraft for a given precipitation type, rate, and temperature. The tables are based on testing of deicing fluids at 1 and 2.5 mm h^{−1} liquid equivalent snowfall rates and at various temperatures. Thus, knowledge of the occurrence of LWE rates less than 1 mm h^{−1} (light), between 1 and 2.5 mm h^{−1} (moderate), and greater than 2.5 mm h^{−1} (heavy) are critical for the proper use of the holdover time tables (Rasmussen et al. 2003). These categories are illustrated in Fig. 1.

The Weather Support to Deicing Decision Making (WSDDM; Rasmussen et al. 2001, 2003) system was developed to nowcast LWE values and is currently deployed at several airports in the United States. WSDDM was also used to provide nowcasts of icing conditions at the Vancouver 2010 Winter Olympic and Paralympic Games (Joe et al. 2010). WSDDM produces nowcasts of LWE values by calibrating the *Z*–*S* relationship, which relates the radar reflectivity factor *Z* to LWE rate *S*, in real time to account for variability in precipitation density. Extrapolated Weather Surveillance Radar-1988 Doppler (WSR-88D) Level-III base reflectivity values are converted to LWE estimates using the *Z*–*S* equation calibrated by previous radar and surface observations. Reflectivity values are extrapolated according to motion vectors estimated by maximizing the cross correlation between subgrids in successive radar data frames.

The Dynamic and Adaptive Radar Tracking of Storms (DARTS) nowcasting method is based on the continuity equation describing the flux and evolution of an observed precipitation field represented by a sequence of radar reflectivity fields. This sequence is modeled as a discrete spatiotemporal linear model that is formulated in the Fourier domain. The motion vector field is then estimated using linear least squares estimation (Ruzanski et al. 2011). This formulation allows for rapid estimates of precipitation pattern motion to be generated.

This paper presents an evaluation of incorporating DARTS into the WSDDM system by replacing the cross-correlation-based radar nowcasting method currently used. Approximately 92 h of data collected using a World Meteorological Organization (WMO) standard reference gauge during four winter precipitation events at the National Center for Atmospheric Research (NCAR) Marshall Field test site and the Front Range Airport (KFTG) WSR-88D near the Denver International Airport (DIA) from 2007 to 2008 were used for evaluation. Results showed an approximate 14% average improvement in nowcasting accuracy, represented by the critical success index (CSI; Wilks 2006) and a corresponding reduction in runtime of two orders of magnitude.

## 2. The WSDDM algorithm

The estimation of liquid equivalent snowfall rate using radar data has been well studied (Langville and Thain 1951; Marshall and Gunn 1952; Carlson and Marshall 1972; Boucher and Weiler 1985; Fujiyoshi et al. 1990). These studies related *Z* (mm^{6} m^{−3}) to *S* (mm h^{−1}) via a power law given by

Several studies noted the variability between reflectivity and LWE rate. Marshall and Gunn (1952) showed a factor of 10 variation in the daily value of the coefficient *a*. Puhakka (1975) observed factor of 5 storm-to-storm variations in the value of *a*. Fujiyoshi et al. (1990) showed that snowfall rate varied by two orders of magnitude for a radar reflectivity of 20 dB*Z* using 1-min radar and snowfall rate data. Rasmussen et al. (2003) showed that both theoretical and empirical *Z*–*S* relationships exhibit a factor of 4 variation in LWE rate for a given reflectivity value.

Super and Holroyd (1998) sought to compensate for variability in the *Z*–*S* relationship in their development of a snowfall accumulation algorithm for the WSR-88D radar. The most significant configurable parameters of their algorithm include the parameters of the *Z*–*S* relationship, the snow–water ratio, and a range correction factor. These parameters can be changed on a storm-by-storm basis.

The WSDDM algorithm accounts for the natural variations in the *Z*–*S* relationship that occur during a storm due to the collective changes in snow crystal type, degree of riming and aggregation, wetness, and size distribution by calibrating the *Z*–*S* equation in real time. The calibration process involves estimating the coefficient, *a*, and exponent, *b*, in Eq. (1) by correlating radar and surface observations whenever new data become available. Nowcasts of radar reflectivity are then converted to nowcasts of LWE using this calibrated equation.

The WSDDM LWE nowcasting procedure is described by the following steps and illustrated in Fig. 2:

Determine the value of the exponent,

*b*, in Eq. (1) according to surface temperature observations,*T*. For pure rain events (*T*> 2°C), the Marshall–Palmer value of*b*= 1.6 is used. For snow events (*T*< −3°C), a value of*b*= 1.75 is used based on the assumption of constant terminal velocity. For mixed events, the value of*b*is interpolated assuming a linear relationship with temperature (Rasmussen et al. 2003).Determine the radar-derived motion vector

**V**over the location of the gauge_{r}*p*by averaging the motion vectors over a chosen area centered over the gauge._{g}Estimate the average fall speed |

**V**| of the hydrometeors according to_{f}*T*. For snow events, a value of |**V**| = 0.9 m s_{f}^{−1}is used. For rain events, a value of |**V**| = 10 m s_{f}^{−1}is assumed. For mixed events, the value of |**V**| is interpolated assuming a linear relationship with temperature (Rasmussen et al. 2003)._{f}Locate the (upwind) location

*p*_{0}of the average reflectivity value*Z*_{0}associated with the LWE value*S*_{0}according to −**V**×_{g}*t*._{f}- Compute the coefficient
*a*in Eq. (1) to make the radar- and gauge-estimated accumulations equal according to where is the LWE accumulation from the snow gauge over the chosen period (30 min was chosen based on results of a preliminary study). Determine the radar-derived motion vector

**V**_{0}corresponding to the point*p*_{0}by averaging the motion vectors over a chosen area centered over*p*_{0}.Locate (upwind) locations of average reflectivity values according to

**V**_{0}and the desired nowcast lead time*τ*by computing −_{r}**V**_{0}×*τ*._{r}Convert the nowcasted average reflectivity value to LWE using the calibrated

*Z*–*S*relationship.Compute the lead time,

*τ*, corresponding to LWE nowcasts at_{g}*p*_{0}by*τ*=_{g}*τ*+_{r}*t*._{f}

A 10 km × 10 km averaging area centered on the point of interest was chosen in all averaging operations based on results of a preliminary study.

Thus, the performance of WSDDM depends on several factors. The quality and resolution of surface temperature, wind velocity, and precipitation data are critical to proper calibration of Eq. (1). Additionally, the accuracy of the radar-derived motion vectors predicated on the resolution and quality of the radar data is important for calibration and nowcasting.

The Tracking of Radar Echoes by Correlation (TREC; Rinehart and Garvey 1978; Rinehart 1981; Tuttle and Foote 1990; Chornoboy et al. 1994; Li et al. 1995) algorithm is the current radar-based nowcasting method used in the operational WSDDM system. To estimate a motion vector, a subgrid is centered over each location in the earlier of the two input grids. The distance between the locations of the maximum cross-correlation value and the subgrid center point is found between the two input grids, and a motion vector is computed (Fig. 3). This process is repeated for every grid point or a desired subset of grid points to estimate the motion vector field.

## 3. The DARTS algorithm

The DARTS algorithm represents the flux of a precipitation pattern represented by a temporal sequence of radar reflectivity fields by (Jacobson 2005)

where *F*(*x*, *y*, *t*) is the sequence of reflectivity fields, *U*(*x*, *y*) is the east–west component of the velocity field, and *V*(*x*, *y*) is the north–south component of the velocity field. Noting the relationship between the discrete Fourier transform (DFT) and the Fourier series expansion (FSE) as well as the derivation property of the DFT (Oppenheim and Schafer 2009), Eq. (5) is represented (Ruzanski et al. 2011) by

where *F*_{DFT}(*k _{x}*,

*k*,

_{y}*k*) represents the 3D DFT coefficients of the discrete observed radar field sequence

_{t}*F*(

*i*,

*j*,

*k*). Additionally, represents the 2D DFT coefficients of the field of estimated east–west motion vector components

*U*(

*i*,

*j*), and represents the 2D DFT coefficients of the field of estimated north–south motion vector components

*V*(

*i*,

*j*). The quantities

*T*and

_{x}*T*, are the lengths of the horizontal and vertical dimensions of the observed gridded reflectivity fields, respectively;

_{y}*T*is the number of reflectivity fields considered for motion estimation; and

_{t}*N*and

_{x}*N*are the maximum harmonic numbers of

_{y}*F*

_{DFT}(

*k*,

_{x}*k*,

_{y}*k*) in the east–west and north–south dimensions, respectively.

_{t}Equation (6) can be implemented in linear form as shown:

where **y** is a vector representing *F*_{DFT}(*k _{x}*,

*k*,

_{y}*k*),

_{t}**x**is a vector representing and , and is constructed such that Eq. (6) is satisfied. The DFT coefficients of the components of the motion vector field can then be estimated using linear least squares estimation by

Further details of the derivation of the DARTS algorithm are given by Ruzanski et al. (2011).

## 4. Data

Archived data collected during four winter precipitation events (approximately 92 h) were considered in this study. Surface data were collected at the NCAR Marshall Field test site, located southeast of Boulder, Colorado (39.9491° latitude, −105.1954° longitude, 1745 m altitude). Precipitation data were collected using a WMO standard reference Geonor T-200B Double Fence Intercomparison Reference (DFIR) gauge. Previous research has shown that without a wind shield, gauges can underestimate precipitation by 50% or more during windy events (Landolt et al. 2004). A study performed by Rasmussen et al. (2001) revealed that the Geonor T-200B in the DFIR shield located at the Marshall site was within 5% of manual measurements. An R. M. Young anemometer was used to collect surface wind velocity data, and a Vaisala WXT-510 was used to collect temperature data. The temporal resolution of all surface gauge data was 1 min. Radar data consisted of Level-III base reflectivity collected by the KFTG WSR-88D radar outside Denver, Colorado (39.7866° latitude, −105.5458° longitude, 1675 m altitude). Data from the lowest elevation angle scan (0.5°) were projected onto a grid with a spatial resolution of 1 km × 1 km covering an area of 200 km × 200 km. The radar was operating in clear-air mode during the duration of the events, yielding a data range of ±28 dB*Z* and an update rate of 10 min. The map of the sensor layout used for this evaluation is shown in Fig. 4.

Table 1 provides details of the events, including statistics of surface observations (i.e., temperature, wind speed, and LWE rate) characterizing each event relative to WSDDM parameters. Surface temperature values and variability relate to the values and variability of the *b* parameter and the hydrometeor fall speed used in the calibration process of the *a* parameter in Eq. (1). The wind speed measured at the gauge site relates to the calibration of the *a* parameter in Eq. (1). The LWE rate statistics and totals give insight into the intensity, variability, and moisture content of the event that relate to the values and variability of the *a* parameter in Eq. (1). Storm duration was determined by the presence of radar echoes and continuous precipitation observations at the gauge site. Collectively, these data quantify a range of conditions under which DARTS was evaluated to yield potential performance improvement of the WSDDM system.

## 5. Performance assessment of the WSDDM system using the DARTS radar–based nowcasting method

### a. Methodology

The performance of WSDDM was evaluated based on the LWE rates categories less than 1 mm h^{−1}, between 1 and 2.5 mm h^{−1}, and greater than 2.5 mm h^{−1}, as illustrated in Fig. 1. The CSI is commonly used to quantify performance in such categorical scenarios and is defined (Wilks 2006) as

where *A* (“hit”) represents the number of correctly nowcasted categories, *B* (“false alarm”) represents the number of overforecasts of the actual rate, and *C* (“miss”) represents the number of underforecasts. Gerapetritis and Pelissier (2004) presented an analysis to characterize the CSI in terms of the false alarm rate (FAR),

and the probability of detection (POD),

The CSI, POD, and FAR range in value from 0 to 1, with 1 being a perfect CSI and POD score and 0 being a perfect FAR score.

The 60-min CSI, POD, and FAR are considered in this study to assess nowcasting accuracy. Whereas the CSI gives insight into the overall performance, the POD stratifies nowcasting accuracy in terms of the ability to detect the correct LWE category. Likewise, the FAR stratifies nowcast accuracy in terms of the aversion (or propensity) toward overforecasting the correct LWE category. The “60-min” prefix refers to the lead time at which the nowcasts are verified. At each time new data is acquired, a nowcast of LWE rate is made 60 min into the future. This nowcast is compared to the corresponding observation made 60 min later. Thus, nowcasts are made up to the last hour of available data.

The central processing unit (CPU) time is used to assess the overall efficiency of the WSDDM algorithm. In contrast to elapsed real time (or wall clock time), which is subject to input/output time delay and all other types of waits incurred by the program, the CPU time is the amount of time used to process the instructions of a computer program. The CPU time is given in seconds and represented (Thimmannagari 2005) as

where IC is the instruction count (dependent on the instruction set architecture and compiler technology), CPI is the average number of cycles per instruction (dependent on the processor architecture), and F is the clock frequency (dependent on the CPU implementation, e.g., nonpipelined, pipelined). WSDDM was run using TREC and DARTS as the radar-based nowcasting methodologies on the same standard Linux compute server.

### b. Results

The CSI, POD, and FAR scores for a forecast lead time of 60 min and CPU times are shown in Fig. 5. Figure 5 shows CSI scores for WSDDM using DARTS versus TREC are higher in three of the four events (11 December 2007, 27 December 2007, and 17 March 2008), the same in one event (4 February 2008), and approximately 14% higher on average. The greatest difference in CSI is seen during the 17 March 2008 event, considered the warmest and wettest event according to the statistics presented in Table 1. For the 27 December 2007 event, the coldest, windiest, and most intense event of the four studied, the differences can be attributed to lower FAR. For the 17 March 2008 event, the warmest and wettest event, the difference can be attributed to a higher POD. Improvements in both FAR and POD explain the improvement in CSI during the 11 December 2007 event, the most moderate event in terms of temperature, wind speed, and intensity of the four studied. WSDDM performed similarly in terms of accuracy using both TREC and DARTS during the drier, less intense 4 February 2008 event.

Figure 6 illustrates the difference in structure between TREC- and DARTS-derived motion vector fields and presents a potential reason for the differences in nowcasting accuracy. Figure 6 shows an example of a reflectivity field and associated motion vector fields estimated using TREC and DARTS methods corresponding to 1227 UTC 27 December 2007. Radar echo motion was observed from the northwest–southeast, and the zero isodop (i.e., the area where the radial is perpendicular to the wind vector) extended from the radar to the southwest and northeast. Quantitative precipitation estimates have been shown to be adversely affected by the zero isodop (Fulton et al. 1998; Nissen et al. 2003), which affects motion vector estimates produced by both TREC and DARTS. The errors in the TREC vector field are apparent as the estimated vectors are not consistent with the observed echo motion and no vectors are present over the gauge location. In this instance, DARTS produces motion vector estimates around the gauge location that are more consistent with the observed storm motion and surface observation. This can be attributed to the interpolation property of the DFT (Oppenheim and Schafer 2009) and the continuity equation upon which the DARTS model is based. Such structure imposes continuity of motion vectors over the entire radar data domain even where no data are available or when data anomalies are present. This provides motion vector estimates over the gauge site used for calibration of the *Z*–*S* equation despite the presence of the zero isodop. Figure 6 thus presents one instance where DARTS-derived motion vectors were clearly better than those estimated by TREC. The occurrence of such differences likely correlates with the difference in nowcasting accuracy scores.

Figure 5 shows that running WSDDM with DARTS reduces the CPU time needed to produce LWE nowcasts by approximately two orders of magnitude for each of the four events considered and overall. These results agree with the theoretical results presented in the appendix and can be attributed to the formulation of each of the two radar-based nowcasting algorithms. DARTS is a Fourier domain–based linear model amenable to computationally efficient algorithms, whereas TREC uses an iterative search method to estimate motion. Figure 5 also shows that TREC exhibits larger differences in CPU time between events (e.g., more than a factor of 2 difference between the 27 December 2007 and 4 February 2008 events). This is because TREC only computes motion vector estimates where data are present and radar data fields were generally sparser during the 4 February 2008 event.

## 6. Summary and conclusions

Nowcasts of LWE values can be used to assist aviation deicing decision-making activities to mitigate costly hazards. The WSDDM system provides LWE nowcasts and is currently deployed at several major airports in the United States. WSDDM produces these nowcasts by calibrating an equation relating radar reflectivity to LWE rate using correlating radar and surface observations in real time. Then, nowcasts of reflectivity are converted to nowcasts of LWE using this calibrated equation. Thus, the radar-based nowcasting method is a key component of WSDDM.

This paper shows the utility of replacing the TREC radar–based nowcasting method with DARTS in the WSDDM system to nowcast LWE rates for aviation deicing decision support. Both TREC and DARTS estimate the motion of precipitation patterns represented by a sequence of radar reflectivity observations. While TREC produces these estimates using an iterative search method, DARTS is formulated in a Fourier domain–based linear model framework amenable to computationally efficient algorithms. Additionally, the structure of DARTS imposes the continuity of motion vectors over the entire radar data domain even where no data are available or when data anomalies are present.

LWE estimates using DARTS were shown to be approximately 14% more accurate on average than those made using TREC in terms of CSI, considering approximately 92 h of data collected at the NCAR Marshall Field test site and the KFTG WSR-88D. The improvement was shown in both POD and FAR scores for three of the four events studied. CPU time was reduced significantly by about two orders of magnitude, which is in agreement with the theoretical results.

These results illustrate the potential utility of replacing TREC with DARTS in an operational WSDDM system. Given the proximity to sensors used in this evaluation (Fig. 4), the results suggest the potential inclusion of DARTS in the operational WSDDM system DIA. WSDDM is currently deployed at DIA, among several other airports in the United States, and DARTS has provided favorable operational emergency decision-making support during severe weather events in central Oklahoma (Ruzanski et al. 2011). Further evaluation using a larger and more diverse dataset should be performed. Such an evaluation should be performed in an operational setting at the desired deployment location to account for the local weather regime, data characteristics, and requirements of operational use.

## Acknowledgments

The authors acknowledge Dr. Roy Rasmussen of the National Center for Atmospheric Research for his strong encouragement to pursue this concept and for his support by providing the appropriate data with which this concept was evaluated.

### APPENDIX

#### Computational Complexities of the TREC and DARTS Radar–Based Nowcasting Methods

The computational complexities of the DARTS and TREC algorithms can be quantified more formally using the Bachman–Landau O notation (de Bruijn 1981). The Bachman–Landau O notation (also called the “Big O” notation) describes the limiting behavior of a function when the argument tends toward a particular value or infinity. The Big O notation thus characterizes functions according to their growth rates. This notation is frequently used to describe an algorithm’s usage of computational resources (Knuth 1997; Corman et al. 2009).

Consider a sequence of *L* data fields of size *N* points × *N* points, where *L* < *N* is on the order of *N*. DARTS performs a forward and inverse 3D fast Fourier transform (FFT) operation, each with complexity *O*(*N* log_{2}*N*) (Cooley and Tukey 1965; Johnson and Frigo 2007). Lower-order computations are performed to construct and retrieve the linear model. A matrix pseudoinverse operation is needed to solve the linear system. Let *M* < *N* be the dimension of the truncated field of FFT coefficients determined by the choice of parameters *N _{x}* and

*N*in Eq. (6). A matrix pseudoinverse operation [using the modified Golub–Reinsch algorithm (Chan 1982)] performed on a single

_{y}*N*×

*N*data field has complexity

*O*(

*M*

^{3}). Considering there are on the order of

*N*data fields used for motion estimation, the complexity of DARTS is given as

*O*(

*N*

^{4}).

Now consider two *N* × *N* data fields between which TREC iteratively searches for points of maximum cross correlation. The complexity of TREC is given as *O*(*q*^{2}*n*^{2}*N*^{2}) (Tsai and Lin 2003), where *q* and *n* (where *q* and *n* are both less than *N*) are the dimensions of the square subtemplate and search area used for pattern matching at each data point, respectively. Both *q* and *n* are determined by parameter selection and are on the order of *N*. The computational complexity of the TREC algorithm is thus given as *O*(*N*^{6}), or about two orders of magnitude greater than the complexity of DARTS.

## REFERENCES

## Footnotes

Current affiliation: Vaisala, Inc., Louisville, Colorado.