Chree Method of Analysis: A Critique of Its Application to Forbush Events Selection Criteria and Timing

The solar–climate investigation has been ongoing for over 200 yr. Definitive conclusions should have, expectedly, been drawn on the subject. This, however, is not the case as is evident in the recent review of the GCR–climate/cloud relationship by Laken et al. (2012). Although the proposed solar–terrestrial links have not been generally accepted, there are numerous factors (e.g., background variability of weather, physical mechanisms, statistical significance tests, data handling, and data quality/selection/smoothing/uncertainties/high-level noise) that could interfere with the detectability of the phenomena, supposing that they exist. While some publications (see, e.g., Haurwitz & Brier 1981; Laken & Calogovic 2013) are dedicated to the implications of statistics of solar–terrestrial analysis, Pittocks (1978) attempted a broad review of a range of limitations associated with solar–terrestrial studies. Data selection is one of the questions pertinent to the present analysis. Pittocks (1978) questioned the validity of statistical significance level associated with representative samples usually selected based on a certain area, time, or specified variables. The problem of selectivity is common among the solar–terrestrial physics community and arises when an investigator’s selection process is directed toward a particular result. Pittocks (1978) has a long list of articles that either carefully or mistakenly selected data biased toward their preconceived ideas. Selectivity and statistical significance are related, especially in cases of marginal statistical significance. Although most authors attempt a justification for their data selection process, Pittocks (1978, and references therein) argued that some of the reasons are apparently unjustifiable. The major characteristic of a biased selectivity in time is unreproducibility of results by other users of the same data. Despite the traditional lopsided appeal to the statistical technique pervading literature on solar–climate influences, Pittocks (1978) opined that data quality and event selection should also be given adequate attention.

Data smoothing and autocorrelations could also have significant influence on the results obtained in solar-weather analysis. Smoothing can be achieved by the methods of normalization, running means, or data filtering. Smoothing usually introduces additional autocorrelation into the data and ultimately reduces the number of degrees of freedom and the attendant greater variability in the correlation coefficients between the data series (Pittocks 1978; Laken & Calogovic 2013). A number of the correlation coefficients, especially the larger ones, could be spurious. In an attempt to illustrate the effects of autocorrelation on GCR–cloud related studies, Laken & Calogovic (2013) compared the confidence intervals computed for individual time points in epoch analysis (hereafter referred to as Method A) with those calculated using the normalization period (hereafter referred to as Method B). Figure 9 of Laken & Calogovic (2013) shows that there is a significant difference between the confidence interval computed based on normalization period (i.e., Method B) alone (and extended over the whole composite time) and those calculated separately for each time of the composite period (Method A). The confidence intervals of Method B are much smaller compared with those of Method A. While articles searching for solar–terrestrial relationships are replete with Method B, Method A (though introduced five decades ago by Schuurmans & Oort 1969) has been attempted by a few investigators (see, e.g., Scott et al. 2014; Okike & Umahi 2019a, hereafter Paper III). Since autocorrelation is one of the common characteristics of geophysical as well as solar data (and usually increased if filtering is introduced), Laken & Calogovic (2013) argued that some of the works claiming significant GCR–cloud correlations have not accounted for quasi-persistence and thus concluded that their submissions are questionable.

Laken & Calogovic (2013) also criticized the small number of FD events usually selected for solar–terrestrial analysis. The Monte Carlo trial methods, commonly used to test the statistical significance of analysis results, are based on the assumption that each data point of a time series has an equal chance of being included in the analysis. However, this is not usually the case, especially in superposition analyses using FD events where the number of Forbush events selected (subsamples) is frequently reduced by timing and various selection criteria used by different authors. The sample size is equally connected to the problems of signal detection and those of signal-to-noise ratio. Harrison & Ambaum (2010) used a large sample of FDs including strong (≤−20%) and small events (≤−3%) to test the impact of FDs on clouds at Shetland. They illustrated an interesting trade-off between sample size/magnitude and noise ratio. When 23 big events (≤−10%) were used, the signal-to-noise ratio was 3.1, whereas it was 2 when only the three strongest events (<−20%) were used. When they employed all 137 events, the signal-to-noise ratio increased to 2.1.

Common among CR data users is the tendency to presume 100% accuracy of data presented on the initial CD-ROM or online via the World Data Center System. Nevertheless, the excellent review of 50 yr CR data by Shea & Smart (2000) shows that a plethora of data handling errors could plague CR data available at various websites. A number of factors, ranging from typo errors, variability of barometric pressure coefficients, normalization factors, inconsistent data values arising from malfunction of a single NM tube, scaling factors, and sensitivity of NM/cutoff rigidity (Moraal et al. 2000), among a host of other unintentional mistakes, could lead to erroneous data. Though these problems have been greatly minimized since the advent of the computer, Shea & Smart (2000) nonetheless point to the need to scrutinize CR data generated by computers.

Definition of selection criteria has been the traditional practice, from inception, among researchers conducting FD-based composition analysis (Marcz 1997). A review of the existing literature suggests that irrespective of the ambiguity of the conditions, any criterion prescribed by a researcher goes without query. This is disturbing, as it might be part of the inherent misunderstanding among investigators in the field. Although we will leave the details for future work, we wish to make a few comments on some of the criteria defined in some publications.

One of the conditions requires the exclusion of FDs accompanied by solar energetic events (Pudovkin & Veretenenko 1995; Todd & Kniveton 2001; Kristjansson et al. 2008; Dragic et al. 2011), though a greater number (e.g., Marcz 1997; Svensmark et al. 2012, Paper II) conducted their analysis without this criterion. Although many researchers attempt to state this as a prerequisite for FD event selections, we note that efforts are yet to be made to confirm the validity of the key event selected vis-à-vis the number of ground-level enhancements (GLEs) that happen within the period. Another criterion, most commonly used, is that of FD magnitude, which is usually determined by a normalization baseline (e.g., Todd & Kniveton 2001; Kristjansson et al. 2008; Oh et al. 2008; Svensmark et al. 2009b, 2012; Kane 2010; Dragic et al. 2011; Okike & Umahi 2019b). Unfortunately, there is no global scale to measure FD magnitude. This is attributed to wide variability in the characteristics of NMs at different locations of the world. It has been observed, for example, that the event of 1991 June 13 is assigned a magnitude of 3%, 7%, 17%, and 30% by different scientists (see Paper II and references therein), suggesting that CR intensity variation that looks like a strong FD at one station might be either a very small FD or even a CR diurnal anisotropy at another. This points to another source of serious bias in FD-related epoch research.

Irrespective of the criterion used, the FD sample size remains surprisingly small. Are the FDs quite limited in number as is generally implied by FD-based Chree analysis, or is it that the problem lies with the researchers’ identification method? The answer to the question becomes obvious when one considers that the IZMIRAN group selected 5900 FDs for the period of 50 yr (1957–2006) using their global survey method, whereas Dragic et al. (2011) identified only 184 FDs for a 42 yr period (1954–1995). While small FDs, detectable with a small baseline of about CR(%) ≤ −0.5, are included in the Belov (2008) catalog, FDs, small and strong events inclusive, cannot be efficiently detected manually. Many FDs that happened within the period of investigation are thus rarely accounted for in FD-based superposition or correlation/regression studies. We feel that verification of these criteria is necessary for statistical significance tests of SEA results.

The CR data used here are sourced from http://cr0.IZMIRAN.rssi.ru/. GCR intensity variations are known to exhibit periodicities, cycles, recurrence, nonrecurrence, and other similar phenomena. The type of analytical transformation applied to observational data depends on the nature of the signal. Figure 1(a) shows the flow pattern of GCR radiation measured at Climax Observatory. The horizontal line represents the average intensity. The sinusoidal wave nature of the signal, driven by the 11 yr solar cycle and the diurnal wave (Kudela et al. 2000), is evident. The long tails pointing toward lower counts are reflective of Forbush effects. The spikes representing solar energetic particles (SEPs) have been removed for full display of the variations of interest. Researchers have developed a number of techniques in an attempt to analyze CR data, some applying general mathematical normalization or common filtering techniques, apparently without recourse to the full profile of the CR time functions (Forbush et al. 1982b; Todd & Kniveton 2001; Singh 2006; Kristjansson et al. 2008; Oh et al. 2008; Harrison & Ambaum 2010; Laken et al. 2011, 2012; Svensmark et al. 2012; Laken & Calogovic 2013; Tezari et al. 2016), and others taking cognizance of the periodic, circular, and sinusoidal wave-like motion of GCRs measured at Earth (e.g., Firoz & Kudela 2007). Tezari et al. (2016), for example, normalized CR data using Equation (2) in their study of latitudinal and longitudinal dependence of the CR diurnal anisotropy. We note that Equation (2) is a variety of Equation (1); the results of these formulae have been discussed in light of our method:

Firoz & Kudela (2007) , on a separate approach, represented CR intensity fluctuation at any time t_i using

GCR_I represents GCR intensity, f(t_i) a wave function, a₀, average intensity, a₁ the daily amplitude, ψ the phase value of the signal, ω the angular velocity, and t_i the time factor. Recalling the expression for sinusoidal motion and simple trigonometric identities, f(t_i) can be expanded to Equation (4), or as Equation (5) following Fourier expression:

While employing the straight line equation based on the normalization/filtering method, Laken & Calogovic (2013) observed that the approach is full of uncertainties, as there is no exact criterion to guarantee that the signal of interest, rather than noise, is not reduced or even destroyed in the process. Following the early indications of Bartels (1935) that meteorological, geophysical, or cosmic data that exhibit periodicities, cycles, and other recurrence tendencies should be subjected to harmonic analyses, we first transformed Equation (5) analytically before further processing. It is assumed here that Figure 1(a) is a Fourier series that can be transformed into a number of sine and cosine terms according to Equation (5). The attempted numerical filtering is similar to that of Barouch & Burlaga (1975), as a number of harmonics representing various signals superimposed on the raw data will be separated as different components. One of the terms, the Fundamental Signal, would have a period equal to that of f(t_i), while others will have shortened periods. Suppose that our assumption and transformations are sustained; then, random variations such as FDs, GLEs, persistent fluctuations (e.g., daily variations), and the qausi-persistence waves (e.g., diurnal waves) inherent in GCR data will be accounted for. While Figure 1(b) accounts for the slow-varying portion, Figures 1(a) and (c) represent the high-frequency signals such as FDs and GLEs. Figure 1(c) is further analyzed (see Okike & Umahi 2019b for further details on the methodology).

Zoom In Zoom Out Reset image size

The blue stars in Figure 1(c) represent the dates as well as the FD magnitudes for the whole period (54 yr) that Climax Observatory has data. While it might be easier to manually detect FDs in CR data spanning a period of a few years, FD manual identification might be intractable given the large volume of CR data as presented in Figure 1(a). In order to amplify the clustered stars for easy comparison with the results of other investigators, we split Figure 1(c) into three, presented in Figures 2–4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1.1. Result Validation/Strong FDs

1. Comparison with IZMIRAN catalog: Our program detected a total of 1571 FDs for the period 1957–2006 using daily CR data, while Belov’s group, which produces the IZMIRAN FD catalog (http://spaceweather.IZMIRAN.ru/eng/dbs.html), calculated a total of 6060 FDs within the same period. A baseline of CR(%) ≤ 0.5 is used for the selection of the 1571 FDs, whereas CR(%) ≤ 1 is used for the 969 large FDs selected within the period. The IZMIRAN FDs are based on hourly data and the dates defined by FD onset, whereas the FD dates indicated in Figure 1(c) are defined by the largest depression in CR flux.

When we investigated hourly data, our program detected a total of 2139 FD events with CR(%) ≤ −3.5 baseline. The increased baseline is in line with the submission of Lockwood (1971) that the ratio of the magnitude of the decrease for a few hours could be about 2:1 when compared with the daily average at any cutoff rigidity; our program uses larger baselines for hourly data than when daily averages are analyzed. When we vary the baseline following Okike & Umahi (2019b), baselines of CR(%) ≤ −3 and CR (%) ≤ −2.5 detected 3905 and 6797 FDs, respectively.

Although the group uses a more accurate event selection approach and CR data from the worldwide NMs (Belov et al. 2005), involving CR physical characteristics such as CR density and anisotropy, Okike & Umahi (2019b) observed that their method and FD data are not widely used by researchers analyzing solar–terrestrial connections (see Section 3 for various techniques common among FD-related studies). The notable large FD event of 2005 January 19 frequently used in many superposition investigations (e.g., Svensmark et al. 2016), for example, was recorded on 2005 January 18 in their list.

2. Strong FDs: The magnitude of FDs, as well as what determines it, has been a subject of interest among CR researchers (see Belov et al. 2001). Traditionally, the magnitude of an event usually expressed as a percentage is an indicator of the size or strength of the event. This quantity depends on a number of variables (e.g., the type of CR monitor, vertical cutoff rigidity, and the detector’s atmospheric depth) and should have unique values for specific locations on Earth. However, a survey of literature shows that for a certain place on Earth, different publications might calculate different values for the same FD event using the same NM data. The magnitude of the largest FD event of 1991 June 13 was calculated as 30% and 7% by two different authors using the same CR data (see Paper II for details).

In a bid to further the investigation on FD magnitude/strength and its implication on atmospheric ionization, some researchers (see Usoskin et al. 2011; Svensmark et al. 2016) attempt to rank FDs and GLEs. Svensmark et al. (2009b) presented FDs ranked according to their impact on the low-altitude atmosphere. The same approach is adopted in Svensmark et al. (2012, 2016). Although ranking of FDs requires assimilation and analyses of all the NM data (a task we will leave for future work) as demonstrated by the IZMIRAN group, it might be interesting to compare the results of Svensmark et al. (2016) and a few other published FDs with our result.

Figures 2–4 are presented for the purpose of such comparisons. The symbols PU1, PU2, ..., LAKEN1, LAKEN2, ..., TOD1, TOD2, ..., CA1, CA2, ..., A1, A2, SV1, SV2, ..., and S1, S2, ..., indicated in Figures 2–4, respectively, stand for FDs taken from the lists of Pudovkin & Veretenenko (1995), Laken et al. (2011), Todd & Kniveton (2001), Calogovic et al. (2010), Ahluwalia et al. (2009), and Svensmark et al. (2012, 2016). Except for PUs in Figure 2 where orders 1, 2, and 3 represent serial number of FDs in the list, the rest denote the magnitude/strength of the events in the respective articles. We note that comparison of FDs requires that the same NM data as those of the publications referred be used. However, Oh et al. (2008) and Paper II showed that strong FDs are global phenomena observed simultaneously at all points on Earth. Such large events are generally selected by investigators searching for solar–terrestrial linkage. Before taking up some case studies where the same CR data are used for comparisons, the strong FD events from the papers referred above are compared with the results obtained using Climax data. Three events PU2–PU4, corresponding to the events of 1969 March 24, 1969 November 9, and 1970 June 18, are shown in their Table 1. The first event in their list (PU1, 1969 March 17) might be a small and nonsimultaneous event, as it was not observed at the Climax station. Pudovkin & Veretenenko (1995) assigned 4.8%, 4.5%, and 2.9% to PU2–PU4, respectively, while the corresponding magnitudes calculated for the events by our program are 1.7%, 0.99%, and 2.7%. Though the difference in sizes of FDs might be attributed to the differences in the FD selection method and/or the differences in cutoff rigidity between the two CR stations (Apatity and Climax), it can be inferred that these three events are strong and hence simultaneous at the two NMs. The number of FDs ≤0.99% between 1969 and 1986 is 519. Except for another selection criterion (exclusion of FDs occurring within the time of solar proton events) imposed by Pudovkin & Veretenenko (1995), the implication of this result is that they analyzed about 13% of the strong FDs that could be investigated within the same period using Climax data. The relative strengths of FDs for the period 1957–1970, as well as the total FDs and the number of strong FDs, are shown in Figure 2. In decreasing order of strength, the four largest FDs happened at 1959 July 18, 1959 July 16, 1959 January 22, and 1960 November 14.

Zoom In Zoom Out Reset image size

Laken et al. (2011) used a comparatively larger number of strong FDs (123). Some of the events are indicated in Figure 3 for illustrative purposes. Judging by the magnitude (−1.07%) of the least event (LAKEN2) in the diagram, the number of events within their period (1978–2006) of investigation is 375. The event of 1982 July 14 is the largest in their least. This also confirmed in Figure 3. Two other large events (1972 August 5 and 1978 May 2) are also conspicuous in the diagram. We note that they did not identify their FD onset date from any particular NM data. Their events dates were rather selected from various sources (e.g., published articles, online website).

Some of the events identified by Todd & Kniveton (2001) are marked in Figure 3. Although they claimed to use large FDs, Figure 3 shows that some of the events they used might be quite small. TOD4, for example, is the event of 1986 November 3. They assigned a magnitude of 4.1% to the event, whereas our program calculated the magnitude of the same FD as 0.55%. TOD2 is also a small event (appearing above the horizontal line in the diagram).

A number of strong FDs reported in the other four papers are indicated in Figure 4. While the magnitudes of AH and CA can be interpreted in light of Figures 1 and 2, the events represented by S and SV are selected by Svensmark et al. (2012, 2016) using a system of FD ranking. Some of these events are displayed here to test the agreement between their method of FD ranking and our approach to calculation of FD magnitudes. SV1–SV5 suggest a good agreement between their ranking method and our event catalog. The same applies to the events represented by S. Except for the event of 1991 June 13, which is strongest in Figure 4 but ranked second to that of 2003 October 31 by Svensmark et al. (2016), the order of strength of S11 and S19 is in line with their ranking. The magnitude of the event on 2001 November 25 (S19) is 1.82%, whereas Svensmark et al. (2016) assign it a magnitude of 39%. The number of FDs ≥ 1.82% between 1987 and 2006 is 122, implying that they investigated about 21% of the strong events that happened within the period.

6.2.1. Pudovkin & Veretenenko (1995)

Although the illustration in the preceding section might have reasonably highlighted the existing gap in the FD key event data selection criteria, timing, and ranking when conducting a Chree analysis, one may argue that some of the glaring differences in the number of FDs selected by the literature referred above and the number selected by our programs could rather be attributed to the type of NM data used. The argument is sound since the characteristics of CR detectors vary widely, coupled with the nonuniform distribution of GCRs over Earth, as well as varying cutoff rigidity (Papers II and III). Pudovkin & Veretenenko (1995), for example, analyzed CR data from the Apatity CR station, whereas Climax data are used in the result presented in Figure 2. Figure 5(a) shows data from the Apatity CR station, while Figure 5(b) shows one of the harmonics of interest. The dates/magnitudes of the three conspicuous events ranked according to their strength from the strongest to the weakest are (1978 May 2)/FD1 (134%), (1982 July 14)/FD2 (132%), and (1972 August 5)/FD3 (130%). It can be observed that these three events are also reflected in Figure 5(b), confirming the sinusoidal fidelity of the present technique. Inspection of their FD list in their appendix shows that these three main events are not included. The increasing order of the strength of the event they selected is denoted with the symbols PU1, PU2, PU3, ..., PU65, with PU65 standing for the weakest event included in the appendix referred.

Zoom In Zoom Out Reset image size

Beyond the three largest events, PU1 is the fourth-largest FD as measured by the red horizontal line. There are, however, many other events with similar magnitudes to PU1 that were not identified by their method. PU2 is much smaller than many other FDs that are not selected. The order of their FD magnitude also does not agree with ours, as it is evident that PU15 is much stronger than PU3. While our code tracked most of their events, a number of their events (e.g., 1973 July 26, 1973 January 9, 1973 September 22, 1973 October 31, 1976 March 30) are not reflected in Figure 5(b). A close inspection of their list reveals that these are among the smallest events they identified. As hinted at earlier, attempts to select low-magnitude FDs manually would be biased by CR diurnal anisotropy if not first isolated as illustrated by Paper III. Out of the 422 large FDs selected by our program for the years 1969–1986, Pudovkin & Veretenenko (1995) identified 65 (about 15%) events. However, they noted that Forbush events superposed in the first 3 days by solar proton events are purposefully excluded in their list. But McCracken (2007) indicated that SEPs occur rarely, so as to justify the exclusion of 85% of the FDs within the same period. Asvestari et al. (2017) had a list of the 16 GLEs recorded by the NM for the 18 yr period. This translates to 0.89 yr⁻¹ and 4% of the large FDs within the period. For the years 1958–2005, they estimated the annual frequency of occurrence of such events as 0.06 yr⁻¹, though Todd & Kniveton (2001) indicated that there are years when it might occur once or even as high as 8 yr⁻¹ in some cases.

Figure 5 (Apatity) indicates that a total of 422 large FDs are detected between 1969 and 1986, whereas 519 strong events are reported when Climax data were analyzed for the same period (Figure 1(a)). In addition to the difference in the number of years used for the two figures, the differences might be attributed to cutoff rigidity between the two stations. In order to illustrate the implication of choice of normalization period on the number of FDs selected within a given time, we included more years in Apatity data. When data between 1961 and 2016 were used, a total of 555 FDs were picked for the same period 1969–1986. The mean CR count for the longer period 1961–2016 (7431.741) is larger than that of 1969–1986 (7431.031), and since our program normalizes CR to the mean value of the period used, a greater number of FDs are selected between 1969 and 1986 when a longer data time series is used. This could also explain the fewer number of events detected by researchers using a few-days running mean as a reference point (see Okike & Umahi 2019b, for details of other factors that determine the number of FDs selected with given CR data).

6.2.2. Laken et al. (2011)

The solar–terrestrial analysis of Laken et al. is both unique and controversial. The number of FDs (123; see their Figure 1) used in the study is much larger than those of previous studies. This should, in principle, solve the problems of small sample size and low signal-to-noise ratio that were earlier neglected. In spite of their tremendous efforts to increase FD sample size, they could not confirm the CR–cloud link. Rather, anomalous cloud changes over the Antarctic plateau, believed to be induced by GCRs, were blamed on solar irradiance. In light of the present work, their methodological approach recommends itself for a reassessment.

Their events were selected from different sources, including literature and websites, for the period 1978–2006. Selecting up to 123 FD key times manually is an onerous task, considering the high variability and hidden periodicities in the CR flux. Although FD key dates have previously been accepted without validation, Okike & Umahi (2019b) recently noticed the need for correct and careful selection and timing of Forbush events. One of the ways of validating key event dates selected for a Chree analysis is demonstrated by Harrison & Ambaum (2010). Following the same method, we present some events in the list of Laken et al. (2011) in Figure 6. The green line with filled circles is the mean. While big FDs can be manually selected without much bias, attempts to select small FDs could be frustrated by the presence of CR anisotropy. A comparison of Figures 6 and 1 of Harrison & Ambaum (2010) would suggest that most of the events in Figure 6 are not FDs. Laken et al. (2011) presented an interesting result of the average of all 123 FDs in their Figure 2(a). Again, a comparison of their Figure 2(a) with some of the individual cases in Figure 6 reveals a bias problem in their technique. If a reasonable number of events in their list are not FDs as is reflected in Figure 6, how does the composite mean in their Figure 2(a) exhibit the sudden depression, maximal depression, and slow recovery characteristics of FDs? We have a similar case in Figure 6. Instead of maximal reduction on the epoch day, some of the events show increases, with some recording maximum depression several days after the key date. These notwithstanding, the mean apparently looks like an FD, though the distinctive features are missing. This is attributable to leveraging, one of the unsolved problems in the Chree method (Okike & Umahi 2019b) of investigation, as we will shortly demonstrate.

Zoom In Zoom Out Reset image size

Figure 7 is presented to show the influence of a single large anomaly on the results of a Chree analysis. It is easy to see that four of the events are not only real but also the strongest FDs in their list, while four events, marked with dashed lines, are probably not FDs, as some of them have their minimum 10 days after the epoch day. Since the SEA method is vulnerable to leveraging that can result from the influence of one or more large anomalies, the mean signal in Figure 7 is weighted toward the large FD values. Okike & Umahi (2019b) suggested that the problem of leveraging can be solved by increasing the sample size or by including the normalization step in superposition analysis, while Laken et al. (2012) concluded that some normalization procedures may result in false positives in an SEA as illustrated in Figures 6 and 7.

Zoom In Zoom Out Reset image size

Though their sample size was significantly increased compared to other reports that used as few as six FDs, 123 FDs, selected for their investigation within a 29 yr period, might look very small. The result presented in Figure 1(b) covers the same duration. The number of large Forbush events selected by our program for the years 1978–2006 is 522, implying that only 23% of the strong FDs that occurred within the time were considered in their analysis. We observe, however, that these numbers will vary between CR stations.

6.2.3. Kristjansson et al. (2008)

This work is a little different from those considered earlier. They selected 22 FDs from the Climax station for the years 2000–2005 and went ahead to validate their results using FDs selected from two other stations, Oulu and Moscow NMs. Events observed simultaneously by the three stations should be strong and global (Oh et al. 2008, Paper II). But selecting FDs visually from three stations is subject to bias, among other limitations. Due to nonuniform distribution of GCRs over Earth, the number of FDs seen by one station may differ significantly from those of other stations. This is because CR intensity variability manifests both longitudinal and latitudinal dependence (Tezari et al. 2016; Okike & Umahi 2019b). FD detection, especially small FDs, requires averaging over a number of stations (Barouch & Burlaga 1975), rather than a filtering method. Assimilating data from a number of CR stations requires empirical orthogonal function analysis to extract the leading principal components (see Paper II for details of the methodology). Figure 8(a) represents the PC1 signal resulting from three CR stations (Climax, Oulu, and Moscow). The correlation coefficient between PC1 loadings and PC1 scores is ≈99 for each of the stations. The high correlation coefficients suggest that PC1 accounts for all the intensity variations in the three stations. The proportion of variance due to PC1 is ≈98%, while PC2 and PC3 account for the remaining 2%.

Zoom In Zoom Out Reset image size

Figure 8(b) is the high-frequency signal containing the high-magnitude FDs. The number of large events detected using a baseline of CR(%) ≤ −1 is 100. Although this numerical averaging over a number of stations is more involved (Barouch & Burlaga 1975, Paper II) than filtering over a single station, it is a unique tool for isolating strong/global FDs than the usual FD magnitude, which varies across Earth. All the FDs as reflected in Figure 8 are strong and simultaneous at the three stations, while the local variations that are station dependent, such as weak or nonsimultaneous FDs, are averaged out (Haurwitz & Brier 1981; Okike & Umahi 2019b). In light of the presented PCA analysis, simultaneity of FD event is a stronger indicator of the strength of Forbush events than the percentage decreases, which vary among authors and between places on Earth. Detailed explanation of the usefulness of this approach, especially with regard to weak and strong Forbush event discrimination in a Chree analysis, is left for future work. We, however, note in brief that Kristjansson et al. (2008) admitted that only 6 out of the 22 FDs selected for their analysis were strong, whereas Figure 8 shows that a larger number of big FDs are observed between 2000 and 2005. Indeed, the period 2000–2005 falls within the time of high solar activity as can be inferred from Figure 8(b). The 22 FDs investigated by Kristjansson et al. (2008) imply using only 22%/15% of the large/total FDs in comparison with the number selected by our algorithms. The authors agreed that their failure to see any clear FD-low cloud connection could be due to the small FD sample employed in their analysis. A larger number of key event times is thus required to confirm their result.

A review of Schuurmans & Oort (1969) and Pittocks (1978) indicates that very many scientists, especially the meteorologists, have been highly critical of the empirical results on solar–terrestrial relationships. Statistical significance tests of epoch results, data selection/quality, data smoothing, and autocorrelation issues are some of the initial sources of doubt on solar–terrestrial analysis. Schuurmans & Oort (1969) responded early to the questions regarding suitability of statistical tests. Haurwitz & Brier (1981) and Forbush et al. (1982a) adopted a different approach to demonstrate the nonapplicability of the traditional t- F-, or p-values in the significant test of epoch results. Several years after its introduction, there seems to be a renewed interest (see Laken & Calogovic 2013; Scott et al. 2014; Paper III) in the statistical method of Schuurmans & Oort (1969), which accounts for the presence of autocorrelation in geophysical data. The detailed illustration of Laken & Calogovic (2013) not only redirected researchers to the correct method of statistical significant test but also tends to lay the age-long problems of incorrect test of significance attendant on solar–terrestrial connection to rest.

However, data selection and smoothing/filtering are yet to be fully considered. This is partly because many solar–terrestrial investigators presume that data selection, for instance, is a preliminary stage that does not require much attention. Another reason, as observed by Shea & Smart (2000), could be due to the fact that some researchers investigating solar–terrestrial links, for example, are not familiar with solar physics and CR topics. Nevertheless, the highly critical work of Pittocks (1978) confirms that adequate attention should be given to event selection. The publication has a long list of authors that dubiously selected data to suit their preconceived results.

In light of these criticisms, the present work attempts to call the attention of investigators conducting SEA studies with the FDs and employing regression analysis to test various solar–terrestrial connections to the need for a detailed and careful method of FD event selection. Ramirez et al. (2013) found that FD is the most spectacular variability in the GCR flux and, as such, accurate detection of Forbus events in a large volume of CR data requires a sophisticated method as opposed to the simple (manual) technique documented in numerous studies (Harrison & Ambaum 2010; Laken et al. 2011). Such careful and accurate methods (e.g., ring station method, global survey method) were first developed by the IZMIRAN group (Kudela et al. 2000; Belov et al. 2001, 2005; Belov 2008), and their FDs have been used in several solar–terrestrial studies (see, e.g., Belov et al. 2001; Belov 2008). While they identified 5900 FDs between 1957 and 2006, for example, Harrison & Ambaum (2010) selected 137 FDs even in a longer period (1952–2006) using the manual method of plotting and calculating the amplitude of each event. A comparison of the large number of FDs identified and used by the IZMIRAN group with the small number of FDs selected by other researchers reveals the pitfalls associated with manual FD event selection. The presented harmonic transformation and the accompanying R program tend to lend credence to the abundance of large FDs as demonstrated by the IZMIRAN global survey method (Belov et al. 2001, 2005). While their method uses FD onset time, our program is based on the time of maximal depression. The large event of 2003 October 31 is, for instance, recorded from 2003 October 28 to 30, with the largest magnitude (29.7%) on 2003 October 29, at the IZMIRAN website. Again, a survey of literature shows that the time of maximal CR reduction, rather than the onset time, is more frequently used in composition analysis (Laken et al. 2011; Svensmark et al. 2012), implying that FDs selected by the present algorithm might attract the attention of solar–terrestrial physicists conducting FD-related Chree analysis.

A note of caution is, however, in place here. FD is such a complex phenomenon that no single algorithm or software can completely account for the physics and knowledge behind it. Although substantial efforts are being usually made by the principal investigator that prepares CR data, Shea & Smart (2000) found that CR data are not always completely without errors. It is thus speculated that some of the events picked by our program could be outliers or other erroneous data arising from data handling. In light of this, the present analysis can serve as a stepping stone to accurate FD identification. In particular, events selected by our program could further be validated with solar wind and interplanetary data corresponding to the calculated event time. Additionally, while our program calculates FD event time and magnitude with great precision, it is incapable of describing the nature of the onset, main phase, or recovery phase of a Forbush event. Thus, an investigator interested in separating different types of FDs can use the key time detected by our program and then plot the data so as to visualize whether a particular event has sudden or slow onset, fast or slow recovery phase, and so on.

The current analysis reconfirms the indication of Okike & Umahi (2019b) that there are several factors that could influence the number and size of FDs identified for a given CR data. One such factor, commonly considered by CR scientists, are pre-defined baselines. Using a baseline of ≥3 (Laken et al. 2011) or ≥5 (Kristjansson et al. 2008), researchers hope to select large or strongest FDs and thus assign event sizes. Kristjansson et al. (2008), for instance, calculated 22% reduction for the event of 2003 October 31. However, the automated method uses a smaller baseline of order 0.5 (for smaller FDs) or 1 (strong events). Figure 1(c) also shows that our program assigns a much smaller percentage reduction to FDs, though the differences in the CRA data temporal resolutions per day, hour, minute, or second also play a significant role in the assignment of FD sizes (Lockwood 1971). The event of 2003 October 31 is 8.8% (Figure 4). While the magnitude of the largest event of 1991 June 13 was calculated as −30% and −7%, respectively, by Gurnett & Kurth (1995) and Ahluwalia et al. (2009), Figure 4 suggests that the magnitude is about 9%. Thus, baselines and FD magnitudes are relative quantities that might assume different values for different researchers/algorithms even for the same FD event. Uniqueness can only be assigned globally to date or time of maximum reduction (onset time may not be the same/simultaneous even for strong events) of FDs, provided that they are simultaneous and strong events. This is also why several investigators use date of maximal depression for SEA studies. The situation is different for weak FDs as is also illustrated by the figures in Paper II.

The result presented in Figure 8 is a preliminary test of FD ranking by our algorithm. Three CR stations are employed here. If all the CR NM data are included in the PCA analysis as attempted by Paper II or the global survey method of Belov’s group, a complete FD ranking catalog might be created. The catalog could account for all the strong and simultaneous Forbush events that have occurred since continuous monitoring of CR flux began. Such a list of large FDs might be useful for investigating the age-long controversial solar–terrestrial relationships.