Predictability Associated with High-Latitude Retrograde Waves in the 1979–80 Winter Season

Retrograde long waves in the higher latitudes of the Northern Hemisphere can episodically attain large amplitudes and sustain coherent phase propagation for 2 – 3 weeks. The potential influence of such waves on extended-range weather forecast has been conjectured but not systematically quantified. Using a set of ensemble reforecast data, this study examined the predictability associated with an extraordinary retrograde-wave episode in the 1979 – 80 winter. Quantified by the anomaly correlation of the 500 hPa geopotential height in the 40 – 70°N latitudi nal band, increased week-2 predictability was found within the subperiod with the presence of coherent retrograde waves. Some individual forecasts made within the retrograde-wave event exhibited the behavior of “return of skills”. The results suggest a future investigation into the relation between the elevated level of anomaly correlation in week-2 and detailed dynamics of the retrograde waves.


Introduction
A distinctive type of coherent westward-propagating disturbances in the higher latitudes of the Northern Hemisphere was identified in observation by Branstator (1987), Kushnir (1987), Madden and Speth (1989), and Lanzante (1990). These disturbances attained their maximum amplitude at approximately 60°N. Vertically, they exhibited an equivalent barotropic structure that extended through the troposphere and lower stratosphere. The retrograde disturbances were identified by their footprint in the zonal wavenumber-1 component of the mid-troposphere geopotential height anomaly, whose period ranged from 2 to 3 weeks. The statistics of the retrograde disturbances based on the wavenumber-1 Fourier component were given in Madden and Speth (1989) and recently updated by Raghunathan and Huang (2019).
The aforementioned classical studies already speculated on enhanced predictability associated with largeamplitude retrograde waves. In particular, Branstator (1987) noted a prominent episode of retrograde disturbances in the 1979 -80 winter that (according to an anecdote in Gilchrist 1986) coincides with a known period of elevated predictability in numerical weather prediction. Lau and Nath (1999) examined potential influences of the retrograde waves on weather over North America. Mo (1999) suggested a connection of the retrograde disturbances to the submonthly variability of precipitation over California. Stan and Krishnamurthy (2019) showed the potential of using a retrograde-wave type of mode for the statistical prediction of temperature over North America on the submonthly timescale. The connection of the retrograde disturbances to the cold air outbreak in North America and East Asia was suggested by Lau and Nath (1999) and Takaya and Nakamura (2005a, b). Given that the timescale of the retrograde disturbances falls within the submonthly (1 week to 1 month) window for which atmospheric predictability is not well understood, the influence of the retrograde waves on predictability is of great practical interests. The subject has not been actively pursued despite increasing availability of archived forecast and reforecast datasets.
With this background, this study used ensemble reforecast data to examine the predictability associated with an extraordinary retrograde-wave event in January-March 1980, which was featured prominently in classical observational studies and called "striking" by Branstator (1987). The event was identified by Madden and Speth (1989) as the longest sequence with coherent westward propagation of the wave-1 component in the 8-year (1979 -1987) data that they analyzed. An updated observational analysis by Raghunathan and Huang (2019) reaffirmed that this event remains one of the most prominent in the last four decades. The case study of this extraordinary episode presented here will serve as the basis for a future investigation on multiple episodes using an expanded reforecast dataset.

Reforecast data and bias correction
We used the NOAA Reforecast Version 1 data (Hamill et al. 2004) for the analysis of the forecast for the 1979 -80 winter. The dataset has been updated to Version 2 (Hamill et al. 2013), which covers the period from 1984 to the present. Our use of the Version-1 data was one time only and was justified by its unique coverage of the extraordinary episode of retrograde waves. The analysis focused on the 500 hPa geopotential height. To circumvent issues with diurnal variation, only 00Z data were used. (Hereafter, "daily data" refers to that for 00Z.) The initial condition in the reforecast, essentially the reanalysis, was adopted as the "observation" to compute forecast errors. The ensemble forecast consisted of 15 members, and each was integrated out to 15 days.
The daily reforecasts for 10 winters, 1979/80 -1988/89, were used to produce the daily (00Z) model climatology of the 500 hPa height. [This is comparable to the setting in Madden and Speth (1989), who used the data from 1979/80 to 1986/87 to define climatology.] Each winter is defined as November 16 -March 16. The data were truncated at Triangular 35 (T35) spectral resolution as in the reforecast archive, but note that the original reforecast model has a higher T62 resolution. [The T35 spectral truncation corresponds to retaining all spherical harmonics, Y n m (λ, φ), for total wavenumber n £ 35, where Y n m (λ, φ) = e imλ P n m (φ), with P n m (φ) as the associated Legendre function and (λ, φ) as the longitude and latitude.] For conciseness, we denote the spatial coordinate (λ, φ) as x. A forecast of 500 hPa geopotential height is denoted as Z (x, τ, t, p, k), where τ is the forecast lead time (τ = 1, 3, 5, …, 13, and 15 days are selected for our analysis), t is the day of the season (t = 1, 2, … are November 16, November 17, …) on which the forecast is validated, p indicates year ( p = 1, 2, … are 1979/80, 1980/81, …), and k is the index to identify an ensemble member (k = 0, −1, 1, −2, 2, …, −7, 7; k = 0 is the control run, and k = −N and N are a pair of perturbed runs with equal but opposite-signed initial perturbations.) The 15-member ensemble mean of Z (x, τ, t, p, k) will be denoted as Z (x, τ, t, p). The observation is denoted as Z OBS (x, t, p), where (t, p) are the day and year of the observation.
The daily climatology for observation, Z -OBS (x, t ), is the 10-year average of Z OBS (x, t, p) over p. The daily model climatology for the ensemble mean with lead time τ, Z -(x, τ, t ), is the average of Z (x, τ, t, p) over p. As a quick illustration, Fig. 1a shows the daily climatology of the spherical harmonic Y 3 0 spectral coefficient of Z -OBS (x, t ) and Z -(x, τ, t ) for selected lead time τ for November 16 -March 16. (We choose a spectral coefficient over a grid-point value as the former absorbs the contribution from the global field. Also, a spectral component with m = 0, such as Y 3 0 , is purely real, which simplifies the presentation.) The daily climatology, created froma 10-year average of the reforecast, is further smoothed in time by a cubic least-square interpolation, as shown in Fig. 1b. The procedure is applied to every spherical harmonic coefficient, up to the T35 truncation, of the geopotential height field to produce the smoothed daily climatology for observation and forecast (the latter as a function of forecast lead time), denoted as Z  OBS (x, t ) and Z  (x, τ, t ), respectively. (Note that Z  (x, τ, t ) is already the 15member ensemble mean.) The difference between the model and observed climatology, is the model bias as a function of day-of-the-season and forecast lead time. With bias correction, the cor- Hereafter, the index p is omitted, understanding that only p = 1 (the 1979 -80 winter) is considered for the rest of the paper.

Definition of anomaly
For the analysis of retrograde waves and for quantifying predictability, one needs to define the anomaly of the observed and forecast fields. The observed 500 hPa geopotential height anomaly for the 1979/80 winter is (again, p = 1 is omitted in the notation for brevity) For the forecast, if no bias correction is applied, the anomaly would be Bias correction is equivalent to redefining the anomaly as the departure from the model climatology at the given lead time, i.e., where Z ¢ BC is the bias-corrected anomaly for the forecast.

Prediction of the long-wave component
The distinctive retrograde-wave episode occurred in the second half of the 1979 -80 winter. Using a set of criteria based on the zonal wavenumber-1 component of the 250 hPa geopotential height at 60°N, Madden and Speth (1989) identified this episode as spanning from January 6 to March 22, 1980. An updated analysis by Raghunathan and Huang (2019) found a similar time span of this event using either 250 hPa or 500 hPa geopotential height anomalies. Figure 2a shows the time-longitude Hovmöller diagram of the observed zonal wavenumber-1 component of the 500 hPa geopotential height anomaly at 65°N from Novem ber 16, 197916, , to March 16, 1980. The red vertical bar drawn at the right indicates the active period of retrograde waves based on Madden and Speth (1989). Comparable Hovmöller diagrams for this episode can also be found in Raghunathan and Huang (2019) and Watt-Meyer and Kushner (2015). Madden and Speth (1989) only documented the retrograde-wave events after January 1, 1980. In the updated analysis, Raghunathan and Huang (2019, see their Table 1) found two minor retrograde-wave events in November 1979. Both did not complete one full cycle of phase propagation for the wave-1 component. For simplicity, in the key statistics to be shown in Fig. 5, we lump those minor events into the period of November 16, 1979 -January 5, 1980, and nominally call it the "inactive period", in contrast to the "active period" of January 6 -March 16, 1980, which falls entirely within the long episode with coherent retrograde waves.
Figures 2b -d are similar to Fig. 2a but for the 15member ensemble mean, with bias correction, of the reforecast with 3-, 7-, and 11-day lead times. (Each of the Hovmöller diagrams is produced by weaving many forecasts together.) Over the active period of retrograde waves, the observed phase of the wave-1 component is well preserved in the forecast for up to 7 days and up to 11 days for the period from early February to mid-March of 1980. Over the inactive period, for example through late December of 1979, the predicted phase of the wave-1 component deteriorates much sooner. Ensemble averaging significantly improves the forecast. The prediction from a single ensemble member is much noisier (not shown), even over the active period. Figure 3 shows an example of the "phase dial" of the wave-1 geopotential height anomaly at 65°N from February 14 to February 29, 1980 (days 91 -106 in the Hovmöller diagrams in Fig. 2). Writing the wave-1 component as Z 1 (λ, t ) = A (t ) sin λ + B (t ) cos λ, its amplitude and phase are defined as R = (A 2 + B 2 ) 1/2 and Λ = tan −1 (B/A), with an appropriate choice of the principal value for the latter. The daily (00Z) values of (R, Λ) are plotted in polar coordinate with R as the radial coordinate and Λ as the azimuthal coordinate, defined as increasing counterclockwise. (A counterclockwise evolution of the phase corresponds to westward propagation.) The black dial in Fig. 3 is the observation, and red and blue are the (bias-corrected) ensemble means of forecast with 3-and 5-day lead times. (Each of the red and blue dials is produced by weaving many forecasts together. It should not be confused as a single long forecast.) This demonstrates the preservation of the phase information in the forecast.
To illustrate the effect of ensemble averaging, Fig. 4 uses ellipses to represent the spread of 15 members of forecast for selected lead times. For clarity, we duplicate the observation (black dial) from Fig. 3. The superimposed red filled circles and ellipses are the forecasts made from February 14 and validated at February 15, 17, 19, 25, and 27, i.e., with lead time τ = 1, 3, 5, 11, and 13 days. The dates are not individually marked on the diagram but can be readily identified by the increase of spread with lead time. The counterparts of the five selected dates for observation are marked by filled circles on the black dial. The center of an ellipse is the forecasted ensemble mean of (R, Λ). The length of the semiaxis in the radial direction indicates one standard deviation (calculated from the 15 ensemble members) in R, and that in the azimuthal direction indicates one standard deviation in Λ. At large lead times, an individual ensemble member can behave wildly and ensemble averaging brings the forecast closer to observation. Figures 5a, b show anomaly correlation as a function of forecast lead time for the 500 hPa geopotential height over the 40 -70°N latitudinal band. The results in Fig. 5a are for the ensemble mean with bias correction and are further averaged over all daily forecasts made through the whole 1979 -80 winter. The black, red, and blue curves are for the full geopotential height field (with all wavenumbers), long-wave component with zonal wavenumbers 1 -2 only, and short-wave components with wavenumbers 5 -9. (Here, we group wavenumbers 1 and 2 together as the wave-2 component shares similar dynamics of long waves to the wave-1 component. In the classical studies mentioned in the Introduction, the two are sometimes analyzed together.) Anomaly correlation is calculated with the area weight (square root of the area of the grid box) multiplied by the geopotential height anomaly. In Fig. 5b, the anomaly correlations for the wavenumber 1 and 2 components (in red and blue, respectively) are shown separately. To show the effect of bias correction and ensemble averaging on the prediction of long waves, in Fig. 5b, we first replicate the red curve (for wavenumbers 1 -2) from Fig. 5a as the solid black curve and then superimpose the shortdashed curve as its counterpart without bias correc- tion. Although bias correction helps in improving the forecast, its effect is relatively minor possibly because the traveling long waves have a small mean after phase-averaging. The long-dashed curve in Fig. 5b is the counterpart of the short-dashed curve but with the former representing the anomaly correlation for single predictions. More precisely, for this curve, we first compute the anomaly correlation for each of the 15 ensemble members before averaging the 15 values of anomaly correlation. This is in contrast to the solid or short-dashed black curve for which the 15 forecast fields are averaged before anomaly correlation is computed for the ensemble-averaged field. One may interpret the difference between the long-dashed and short-dashed curves as the effect of ensemble averaging. The noticeably lower value associated with the long-dashed curve illustrates the significant benefit of ensemble forecast for propagating long waves. For the rest of the paper, all results shown are those from the 15-member ensemble mean with bias correction.

Contrast in predictability between active and
inactive periods Figures 5c, d compare the anomaly correlations averaged over the whole 1979 -80 winter (solid lines), inactive period (November 16, 1979-January 6, 1980, and active period (January 6 -March 16, 1980, long-dashed lines) of retrograde waves. Figure 5c shows those for zonal wavenumber-1 (black) and wavenumbers 1 -2 (red) components, and Fig. 5d all wavenumbers (black) and short waves with wavenumbers 5 -9 (red). In all cases, anomaly correlation is higher over the active period. Interestingly, enhanced predictability of the short waves is also found within the active period of retrograde long waves, although further studies are needed to establish a causal relation between the two.
To further demonstrate the contrast in predictability between the active and inactive periods, in Fig. 6a, the anomaly correlation for zonal wavenumbers 1 -2 of the 500 hPa geopotential height (over the 40 -70°N latitudinal band as before) is shown for the individual sets of daily forecasts with the initial condition step- (Each curve in Fig. 6a is for the ensemble mean of 15 forecasts, with bias correction as usual. Averaging over all curves would recover the red curve in Fig.  5a.) A curve is red if the anomaly correlation exceeds 0.5 on day 11, and blue otherwise. Red curves reside overwhelmingly within the active period, indicated by the horizontal red bar at the top of the panel. Figure 6b shows the daily (00Z) values of the phase of the observed zonal wavenumber-1 component of the 500 hPa geopotential height anomaly at 65°N. Figure 6c shows the amplitudes of the wavenumber-1 (black) and wavenumber 1 -2 (blue) components. The latter is defined as A = (|R 1 | 2 + |R 2 | 2 ) 1/2 , where R 1 and R 2 are the amplitudes of wave-1 and wave-2 components. We note that enhanced predictability is associated with the presence of coherent westward phase propagation, rather than just high amplitudes, of the wave-1 component. Figure 7 is similar to Fig. 6a but for the anomaly correlation of the full 500 hPa geopotential height anomaly with all wavenumbers. With the inclusion of short waves, anomaly correlation is reduced overall, but the contrast in predictability between the active and inactive periods is still evident. We put a red dot below the abscissa to indicate an occasion when the anomaly correlation for a forecast exceeds 0.5 at τ = 11 days. (The dot is placed at the date when the forecast is validated. For brevity, in Fig. 7, we omit the labeling of the abscissa, which is otherwise the same as that in Fig. 6.) The active period is abundantly populated with red dots compared to the inactive period.

Spatial patterns of long waves
Because the most significant increase in predictability is associated with the long waves, we next show some forecast maps (all ensemble-averaged and biascorrected) of the long-wave components of the 500 hPa geopotential height. The forecast in Fig. 8, initialized on February 14, 1980, is randomly chosen from the red cases within the retrograde-wave event in Fig.  6a. The wavenumber 1 -2 component of the geopotential height at the initial time is shown in the top left panel. Below it, the left and right columns compare the observation (left) and forecast (right) for February 17, 21, and 25, corresponding to lead times of τ = 3, 7, and 11 days. For this case, the anomaly correlations from Fig. 6a are (0.966, 0.940, 0.865, 0.794, and 0.780) for τ = (3, 5, 7, 9, 11) days. At τ = 11 days, the value of 0.780 is ranked 14th, or just clears the upper tercile, of the 42 red cases within the retrograde-wave event in Fig. 6a. Figure 8 demonstrates the semblance of the predicted long-wave patterns to observation into week-2, although at τ = 11 days, the amplitude of the long waves in the forecast is weaker than its counterpart in the observation.
We next revisit a case singled out by Branstator (1987) as a folklore example of an episode with unusually high predictability (overall, not restricted to long waves) in operational numerical weather prediction. (Branstator cited Gilchrist (1986), which, however, does not contain the detail of the original forecast.) The forecast, shown in Fig. 9 in the same fashion as Fig. 8, is initialized on January 18, 1980. The anomaly correlations for wavenumbers 1 -2 over 40 -70°N are (0.923, 0.824, 0.708, 0.716, and 0.662) for τ = (3,5,7,9,11) days. The corresponding anomaly correlations for all wavenumbers (related to Fig.  7) are (0.881, 0.762, 0.627, 0.594, and 0.547). The observed pattern of the long waves is well preserved in the forecast into τ = 11 days. Although retrograde Red curve and red dot (placed on day 11) indicate the case with AC > 0.5 on day 11. For brevity, we omit the labeling of the abscissa, which is otherwise the same as that in Fig. 6.
waves generally attain a maximum amplitude in the higher latitudes, in this case, the semblance between forecast and observation at τ = 11 days even extends to the lower latitudes (south of 40°N) over North America. This reminds us previous statistical studies that related the retrograde waves to submonthly variability of precipitation (Mo 1999) and temperature (Stan and Krishnamurthy 2019) in the lower latitudes of North America. Although extended-range forecasts with an unusually high anomaly correlation or "return of skill" are known in folklore and some of which could be attributed to random coincidences (e.g., Anderson and van den Dool 1994), we hope that the results here will revive interests in further investigations into nonrandom relations between enhanced predictability and retrograde waves. Raghunathan and Huang (2019)   noted in their potential vorticity (PV) analysis that some strong retrograde-wave events are associated with a poleward extrusion of low PV air over the North Pacific, followed by wave breaking and vortex shedding, and then a westward migration of the anti cyclonic vortex. Such a sequence of distinctive dynamical features is found in the case in Fig. 9.
Using the ERA-Interim data (Dee et al. 2011), Fig. 10 shows the daily maps of PV on the 315 K isentropic surface over the Asian-North Pacific sector for the four dates in Fig. 9. A detached anticyclonic vortex can be identified in the map for January 29. In future work, it will also be interesting to clarify whether the elevated anomaly correlation is associated with a particular stage of the dynamical evolution as illustrated in Fig. 10.

Concluding remarks
This study provided further evidence of enhanced week-2 predictability in numerical weather prediction during a long episode of coherent retrograde waves.
In addition, we quantitatively demonstrated how prevalently the elevated predictability for daily forecast occurred within the retrograde-wave event. Our conclusion is drawn from the reforecast in only one winter, which nevertheless contains one of the most extraordinary episodes of retrograde waves. Motivated by this case study, a follow-up investigation that will analyze multiple episodes of retrograde waves from 1984 to the present using the NOAA Reanalysis II data is underway. As noted by Raghunathan and Huang (2019), only a few retrograde-wave events after 1984 are comparable to the 1979 -80 event in terms of phase coherence and duration. Moreover, the retrograde-wave events documented in Raghunathan and Huang (2019) have a wide range of period and duration. Our forthcoming work will address the issue of weaving those different retrograde-wave events into a unified framework for analyzing the anomaly correlation.
This paper investigated how predictability is influenced by retrograde long waves while the dynamical mechanisms of the waves are not analyzed further. [A small number of works on the dynamics can be found in Branstator and Held (1995), Huang and Robinson (1995), Lau and Nath (1999), and Polvani et al. (1999), also see a brief survey in Raghunathan and Huang (2019).] Nevertheless, analyzing the products of numerical weather prediction as performed in this paper might provide an alternative way to understanding the dynamics of the retrograde waves from an initialvalue-problem point of view.