2. Materials and Methods
While the work in Lorenz [
1] was limited by a lack of very good APs, in Lorenz [
3], better APs were found by using GPH data from ECMWF’s operational model. The later study’s results confirmed his earlier estimate of a 2.5-day doubling period for the initial rate of error growth between a generalized model and the true atmosphere. The chief goal of the present study is applying a modified version of his method to different low wavenumber Fourier components in order to find a range of system predictability estimates.
The logical connection between our numerical findings and the conclusions we may draw from them about system predictability has been based upon Lorenz’s work, but our approach is not identical to Lorenz’s. Our premise is that a large number of fairly good APs will describe system predictability better than a few exceptionally good APs.
Although we begin by generating a roster of APs in a manner very similar to Lorenz [
1], we have applied several new ideas to Lorenz’s method. We propose that, within the space of all two-state comparisons and how they evolve in time, the predictability of this system is well characterized by following the APPs formed from many fairly similar, but not necessarily exceptionally similar, APs. Our use of 100 m in GPH as the difference threshold for choosing APs is a choice that deserves further explanation.
It was by finding a compromise between two opposing considerations that we arrived at the AP set GPH difference threshold value (TV) criterion of 100 m. The first consideration was that the range of difference values between a given AP set’s average difference value and the Esat value, the ‘error saturation’ value at which similarity is lost, should be as great as possible. This is so that our APP curve, built from the chosen set of APs, may traverse as much of the hypothetical ideal curve’s difference range (which ostensibly would extend to difference values close to zero) as possible, given our data. This consideration would have us choose only a few—possibly only one—of the lowest difference-valued APs.
The second consideration is that any very small set of APs will exhibit a substantial amount of variability in the resultant APP curve of averaged values. The effect of this type of variability is that such an APP curve will not rise smoothly to cross the Esat line at such a point as to give a clear indication where to demarcate the limit of predictability. For example, the ‘TV 75 m’ curve in
Figure 1 (darker blue datapoints) amply shows this kind of excessive variability, to the extent that we must question whether the apparent limiting value resulting from that curve (i.e., 20 days) is an accurate reading. This second consideration therefore encourages us toward choosing an AP set sufficiently large to suppress this variability and to smooth out the APP curve. As can be seen by comparing the different curves in
Figure 1, the increase in AP set size results in smoother curves.
The minimum AP difference value found in the background comparison set is 72.7 m. The Esat value is 155.2 m. What value would be the optimum TV? In attempting to qualitatively reconcile these two considerations for choosing TV, we chose 100 m because this was judged to be the lowest difference value that gave us sufficiently large AP set sizes to smooth out the consequent APP curves. This choice gave us an AP set of size 22,287. As we shall see, due to the additional narrowing choices that our different tests required, this number will be reduced (in some cases to sets of only several hundred APs), thereby further limiting the AP set size from which to form APPs. While not perfectly smooth, we judge these resulting APP curves to be sufficiently smooth such that they yield reasonably sharp predictability limit results. A larger AP set would smooth these curves out even more, but we would thereby lose even more range at the low end by starting at a higher average difference value on day 1. This ‘low-end’ range of the average AP set difference values covered by the different TV choices is exactly as seen in
Figure 1 between the datapoints of day 1.
We admit that the 100 m TV criterion is to some extent an arbitrary choice. However, we judged that small changes in that value—perhaps 98 to 102 m—and the different AP sets chosen would not greatly affect our predictability limit results, but that large changes in the criterion value outside that range would no longer reconcile well the two considerations just discussed and would have deleterious effects on the predictability results. On this basis, we assert that the 100 m choice is close to optimum.
Our use of a large set of APs is a significant departure from Lorenz. We propose that we should not limit ourselves to only pursuing the very ‘best’ APs in assembling our roster. Rather, the loss of similarity as a general principle is driven by the accumulation of differences between many somewhat similar, but not very similar, states. The use of many APs, and the APPs built from them, give us a better view into the system as a whole.
We can construct and interpret the behavior of different ensembles of APPs according to different comparison regimes. These regimes can be narrowly chosen to characterize different, specific aspects of the system. We have chosen to follow the APPs constructed from single low wavenumber Fourier components of the 500 hPa GPH field, over the wavenumber range 2–8. We show by comparing the spatial extents of examples that some of these particular components are closely related to Rex blocks.
By following the behaviors of APPs formed from low wavenumber Fourier components and deriving system predictability limits thereby, we are testing the hypothesis that predictability limit is inversely related to component wavenumber. This hypothesis is a formal statement of the intuitively attractive idea that large, persistent structures, such as Rex blocks, should possess an inherently greater range of predictability in time than smaller and more ephemeral structures. It also allows us to draw a direct line from the reanalysis data to Rex block predictability.
We argue that wavenumber components 2–5 are a suitable representation for the large structures—that is, the strong ridging of Rex blocks—with which we are primarily concerned.
Figure 2 left panel shows both the 40° N 500 hPa GPH signal and two other curves for the strong eastern Pacific ridge of 23 January 2014, one formed from the sum of Fourier component wavenumbers 1–8 and one from wavenumbers 2–5. This ridge is seen in the left panel as the prominent peak centered on −130° longitude (also refer to the left panel of Figure 4 for a broader illustration of this ridge). The wavenumber 1–8 curve is quite close to the original GPH curve, and the wavenumber 2–5 curve, while clearly missing some higher-resolution details, still covers the gross structure quite well.
Panel (b) of
Figure 2 shows the difference between the GPH signal in panel (a) and the sum of wavenumbers 1–8 (blue curve) and then the difference between the GPH signal and the sum of wavenumbers 2–5 (red curve). While the maximum difference seen in the GPH minus wavenumbers 2–5 curve is somewhat large, overall, the differences are relatively small, when compared to the amplitudes seen in
Figure 2 (i.e., <13%). The mean difference values for these curves are 19 m for GPH minus wavenumbers 1–8 and 85 m for GPH minus wavenumbers 2–5. The 85 m mean difference value is only a small fraction (<13%) of the maximum range exhibited by the GPH curve of 670 m. For the purpose at hand, resolving big ridging episodes, we judge that wavenumbers 2–5 are the key wavenumber components.
Another way of illustrating this issue is to compare the longitudinal extent of one peak each of wavenumbers 2–5 with the ridge shown by the GPH signal centered at longitude −130°.
Figure 3 shows this comparison, wherein we see that the longitudinal extents of wavenumbers 3, 4, and 5 come closest to the extent of the ridge. The blue GPH signal spans 50 degrees of longitude along the red dashed mean value line, from −100° to −150°. Wavenumber 2 spans 90 degrees, wavenumber 3 spans 60 degrees, wavenumber 4, 45 degrees, and wavenumber 5, 36 degrees. Wavenumber 2 is clearly much broader than this ridge, but the extents of the other wavenumbers 3–5 nicely bracket the size of the ridge.
The data used were obtained from the NOAA ESRL-PSD ‘NCEP-DOE Reanalysis 2’ data set, daily mean values [
17]. These data are somewhat coarse in their spatial increments of 2.5° in both latitudinal and longitudinal directions, giving 144 data points around any given latitude circle and 37 data points along a given longitudinal meridian from the equator to the North Pole, inclusive. We used 90-day winters (DJF; 29 February excluded) for the 38 winters 1979–1980 through 2016–2017. Our study focused on these data’s portrayal of the 500 hPA geopotential heights field, but one part of it (generating a roster of Lorenzian analogue pairs, as described below) used the 850 and 200 hPa fields as well.
Figure 4 shows two examples of the 500 hPa geopotential heights over the area from −170° to −110° western longitude (extending roughly from Hawaii to Utah) and from 20° to 80° north latitude. Note in each panel the presence of a prominent high-pressure ridge over the eastern Pacific Ocean. As this geographic area figures prominently in our study, we will refer to it as area A1. The choice of these days was due to the particular similarity (based on Lorenzian RMS differencing, evaluated over the northern hemisphere) between the two GPH data fields of 23 January 2014 and 9 December 1986. Although this pair of day-states showed a marked visual similarity in the plotted fields, for the three days on either side, we observed only vague visual similarities at best. Such is the rare and ephemeral nature of close RMS similarity in blocking structures: it is quickly lost, as Lorenz noted in his studies of how similarity diminishes between analogue pairs [
1,
3].
As we proceed with describing in detail how we used these data, we wish first to make clear how our process involved several steps that progressively narrowed the number of the two-state comparisons used. The fundamental divide among these steps was between (1) our establishing a large set of APs, in which we adhered closely to Lorenz’ RMS comparison protocol, and (2) our introduction of Fourier decomposition to act upon this set of APs, in order to measure how system predictability varied among combinations of different wavenumber Fourier components and the presence or absence of ridging in A1.
The differencing method used by Lorenz [
1] (p. 638) was used by us to define the extent of similarity, or difference, between the states of any two chosen day-states. The reanalysis GPH data at three pressure levels of 850, 500, and 200 hPa were used, over the northern hemisphere, the equator included. Point-by-point RMS differencing over the chosen data grid was averaged for each pressure level, giving one difference number, in units of meters. We used the following process to calculate the RMS difference between two day-states. A value for
Dijp, the RMS difference for one pressure level, was given by
where
i and
j were two day/states,
p is the pressure level,
Zipk is the GPH value for day
i, gridpoint
k, and the summation was taken over
S gridpoints. The three numbers
Dijp for the three pressure levels
p were averaged, to give just one representative number for the cumulative difference between the two days being compared. This difference number was calculated for every pair of day-states in our 38-winter domain (pairs within the same winter excluded; see below), and it resulted in 5,694,300 2-day comparison values from which 22,287 comparisons were found to be analogue pairs (the AP100 set) by their difference values falling at or below the 100 m GPH comparison TV discussed above. The AP100 set is therefore only defined by the RMS differencing process shown above, consistent with Lorenz’ method. Note that the GPH values used at this step are the pure reanalysis values and are not filtered by Fourier decomposition.
We also chose to disallow intra-winter comparisons—that is, no APs were included in which both dates were chosen from the same winter—so as not to introduce the higher correlations of states due to persistence caused by close temporal association. No additional stipulations, such as the presence or absence of ridging over the eastern Pacific, were introduced at this step. Lorenz used a weighting to adjust his two-state comparisons, to account for drawing days from different seasons, since different seasons show differing variability in GPH range. As we have confined our inquiry to 90-day winters, we have omitted this type of adjustment. We have assumed instead that such variability within the winter season (i.e., between the forcing regimes expressed by 1 December versus 15 January versus 28 February) is not large enough to affect our results. Up to this point, we have not used any Fourier decomposition.
As mentioned above, the AP100 set of analogue pairs, calculated by a purely Lorenzian scheme of RMS differencing, is for the present study merely a starting point. The chief innovation of our study lies within our application of different comparison regimes to the initial AP100 set. The choices of comparison regime were based upon the parameters of wavenumber, phase angle, and amplitude that only become available upon the imposition of Fourier decomposition. By choosing comparison regimes based upon these parameters, we were able to construct APPs whose predictability results illustrated specific qualities of system predictability related to large-scale structures such as Rex blocks.
We now describe our methods for creating APPs that are based upon the 500 hPa GPH data from just one latitude circle at 40° N using Fourier component wavenumbers 2 through 8. Each component was treated separately, to make its own small family of analogue pair progressions.
The choice of using wavenumbers 2–8 was based upon two considerations. The first was that their amplitudes constitute most of the total wave amplitude. For example, on 23 January 2014, the total amplitudes of wavenumbers 2–8 constituted 67% of the total signal amplitude. The second consideration was that the half wavelength extent of these components’ waveforms covered the longitudinal range we expect in regard to the size of ridging in A1. However, as we argued at the beginning of this section, it seems that wavenumbers 2–5 determine the structure of Rex blocks (
Figure 2 and
Figure 3). Wavenumbers 6–8 were included for completeness to show that they do not contribute much to the analysis and thereby bracket the contribution of wavenumbers 2–5 by contrast.
We decomposed the AP100 set as follows. For any day of interest, within our temporal range, the chosen latitude of 40° N resulted in a subset of the 500 hPa GPH reanalysis data with 144 entries. Fourier decomposition was performed on this set. Fourier analysis theory states that, for a given set, there is only one unique component for any single wavenumber (e.g., Ref. [
18]). Therefore, for each wavenumber, the triplet of frequency, amplitude, and phase angle values are unambiguously designated by the Fourier decomposition process.
Wavenumber zero is the set mean value, a constant value around that latitude circle, and thus, it has an ‘amplitude’ defined by the Fourier decomposition process to be the mean value of the GPH signal; its frequency is zero, and its phase angle is undefined. The rest of the Fourier components numbers 1–143 are sine waves with zero mean value, with frequency, amplitude, and phase angle values provided by the decomposition. The variability of that set of 144 values lives within the components 1–143. Or, more properly, we should stipulate that due to the Nyquist frequency of 144/2 = 72, the variability lives in wavenumbers 1 through the ‘fold’ component defined by the Nyquist frequency, wavenumber 72; the higher wavenumber, trans-folds components’ (73–143) amplitudes to exactly mirror their lower-fold (1–71) conjugates. In keeping track of the component amplitudes, one must account for the splitting of amplitudes that this mirroring causes. Since these fold conjugate pairs are identical in amplitude, we must multiply any given wavenumber amplitude from this decomposition process by two to arrive at its correct value, except for wavenumber 72 itself [
18,
19].
In our analysis, each wavenumber family had three members, according to three different initial parsing regimes. These regimes were as follows: S1, parsed for ridging in both amplitude and phase angle (as defined in the paragraph below); S2, the analogue pairs left over from the S1 parsing; and S3, the total set of analogue pairs. Numerically, S1 + S2 = S3, and whereas S1 and S2 vary in size, S3 had a fixed size of 8502 analogue pairs. Note that S3 at size 8502 pairs is much smaller than AP100′s size of 22,287 pairs because we also imposed a separate parsing onto each set to address a temporal issue. We did not want our 21-day progressions extending past February and into March, and thus, we imposed the restriction to only allow analogue pairs that would satisfy this stricture; S1, S2, and S3 all conform to this restriction. We refer to this expedient as ‘date truncation’.
The initial parsing for ridging to form regime S1 used the restriction that, to be part of the set, both members of an analogue pair had to satisfy ridging in A1, and to have an amplitude equal to or greater than the average amplitude of that wavenumber over our 38 winters. Ridging in A1 was defined by choosing phase angle ranges such that a wave peak for that wavenumber would lie within the longitude range of −160° to −120°. This range was chosen as a suitable subset of A1, in that we decided to exclude from our ridging definition those wave peaks lying near (within 10° of) the boundaries of A1. Set S2 therefore was composed of exactly those analogue pairs that failed the S1 parsing definition. Our choice of the S1 parsing was designed to be generally consistent with the structure of the high-pressure ridges of Eastern Pacific Rex blocks, in terms of wave peak position within A1 and wave amplitude.
From the sets of analogue pairs defined by each initial parsing, we built 21-day progressions of the difference values between the analogue pairs as each stepped forward in time. Day 1 shows the difference values between the states of the analogue pairs themselves; day 2 shows the difference values between the states one day past each analogue pair member; and thus, on up to 21 days. Since the analogue pairs were chosen specifically due to their Lorenzian state similarity, we expected that the progressions should show a gradual decrease in similarity in the increase in difference values. Although the resulting curves are not monotonically increasing, they in general follow an increasing trajectory toward the Esat lines. We suspect that larger numbers of analogue pairs would tend to smooth the curves even more and reduce the effects of day-to-day variability.
A progression length of 21 days was chosen due to two considerations. First, we wanted there to be the possibility of trans-Lorenzian predictability limits—more than 14 days—to present themselves. Second, we did not want the progression length to be so large as to drastically limit the analogue pairs chosen. For if we had chosen a progression length of, say, 45 days, but required our progressions to not extend into March, then only those analogue pairs, both of whom started in December and the first half of January, could be chosen. An avenue of future research may be to chart the effects of employing longer progression lengths.
We note at this point that the phase angle and amplitude differences we are using to construct APPs are a natural evolution of the concept of error as used by Lorenz in his 1969 study. In that work, the pattern of difference growth between an analogue pair is analogous to the error that grows between forecasts and the true atmosphere. There is no reason that measuring the accumulation of this error should be confined only to RMS differences of the GPH field. In our work, we have constructed different distributions of two-state difference values for the chosen fields among the many day/states of our temporal domain. The subsequent patterns of variation within these distributions provide different views into the predictability of the system as a whole.
An extensive justification of this strategy seems unnecessary, since it is well known that many different fields are calculated separately for standard forecast model outputs. That is, within the predictive realm of atmospheric science, multiple kinds of inquiries and results are embraced as a matter of course. Our efforts should be seen in this same general light, of choosing to use specific analytical tools of Fourier analysis to focus narrowly on differing aspects of the atmospheric system.
We now turn to an application of the same method, but instead of just one latitude circle, we expanded the analysis to include a broad swath of the northern hemisphere, from 20° to 80° N. As in the previous section, our analysis included the wavenumber range of 2–8.
Our second APP study centered primarily on phase angles. The importance of these phase angles can best be understood when we consider which aspect of waves it is that most characterizes an Eastern Pacific ridge: its longitudinal position, which when decomposed, is purely a function of the component phase angle. The method devised to include 25 latitude circles did not lend itself to including amplitude information, and our attempts proved problematic. We hope a future iteration of this work can include both phase angle and amplitude information, when expanding the analysis from single to multiple latitudes. It would for example be an improvement in the experimental design to impose an improved initial parsing regime to exclude those analogue pairs that have primary wavenumber components of low amplitude, although lying in the correct phase angle ridging range.
We chose the latitude swath of the 25 increments from 20° to 80° N. In our data, these 25 complete latitude circles each have 144 longitudinal increments. Each of these latitude circles was decomposed as an independent signal containing 144 values, yielding 144 components (72 below-fold components) by Fourier decomposition. The choice of using 20–80° N as the latitudinal range of the large swath study was arrived at via the following consideration. We wished to include a latitudinal range sufficient to cover the extent of the varying structure of Rex blocks. This must include most of the northern hemisphere, but we did not want the complication of including tropical and equatorial processes, nor those peculiar to the region immediate to the north pole. While this study investigates processes of the mid-latitudes, its design should not exclude potentially important information that may lie off of the latitudinal fringes of a narrow reading of mid-latitude extent—perhaps, 30–60° N. As shown in the two examples of
Figure 4, the structure of Rex blocks can extend over a greater latitudinal range than this. Our inquiry is into the nature of large-scale patterns, and thus, being more inclusive in our choice of latitudes seemed appropriate for this study. How variability and predictability might change with latitude may form the basis of future research.
For the 25 latitude circles, the data for one day then yielded 25 triplets for each of the 144 components. As in the previous section, we developed each progression sequentially for each wavenumber 2–8, and thus, for each wavenumber, the only changing values were those of amplitude and phase angle. In effect, a day’s full decomposition of 25 latitudes by 144 component triplets was reduced to an array of 2-by-25 values: the given wavenumber’s amplitudes and phase angles for the 25 latitude circles. In order to compare day-states of the analogue pairs and their progression, we devised an averaging scheme to find phase angle centroid values relative to the center of ridging in A1, at longitude −140°.
To give a sense of the general pattern of these progressions, an example is shown in
Figure 5, which shows five time slices along an APP example from the first, single latitude circle study. In a gradual ‘dispersal of the blob’, we can see that the differences in the aggregate between all the APs slowly increase with time, as one would expect. The pattern of increase is obvious over the first three slices, but not as obvious thereafter. Within each progression, the two-way mean values of all the difference values from each day are represented by a large black ‘X’. It is the daily sequence over the 21 days of the progression provided by these two-way mean difference values that yields an estimate of system predictability, as shown in
Figure 6 and
Figure 7.
Figure 6 and
Figure 7 show how the estimate of predictability is made: the progression day at which the Esat line is crossed defines the time value at which predictability has been lost, and thus, in
Figure 5, 19 − 1 = 18 days is the predictability estimate. In a number of the progressions, the progression curve failed to cross the Esat line, its values lying below it for all 21 days, as in
Figure 7. For these progressions, we have termed their predictability estimate as ‘21+’, meaning that according to the evidence of our results, the limit of predictability exceeds our stipulated progression time range.
Table 1 shows that these 21+ day estimates chiefly occurred for wavenumbers 2 and 3, but a few also occurred for some entries of wavenumber 5.