On the probability of occurrence of extreme space weather events

1. Introduction

In analogy with terrestrial weather, “space weather” is the state or condition of the space environment surrounding and within the Earth's magnetosphere. And, just as terrestrial events such as tornados and hurricanes can have devastating effects on Earth, severe space weather events can produce a host of consequences that impact society, including the disruption or loss of space-based assets, such as spacecraft [e.g., Reeves et al., 1998], and terrestrial assets, such as electric power transmission networks [e.g., Bolduc, 2002]. As our society becomes progressively more dependent on technology, assessing the impact of, and ultimately predicting space weather events, and particularly the most extreme ones, will become ever more crucial. In this study, we seek to answer, or at least address the question: What is the probability of such an extreme event occurring within the next decade, or 100 years?

“Extreme” means that which is far removed from the mean or average. But, while an event might be extreme with regard to one parameter, it might only be average relative to others. Moreover, there are an almost limitless number of parameters that can be used to describe the severity of a solar event. Beginning at the Sun and moving toward the Earth, these include, but are certainly not limited to: (1) soft and hard X-ray flares; (2) coronal mass ejection (CME) speed; (3) shock strength; (4) Solar Energetic Particle (SEP) fluxes, including Solar Proton Events (SPEs); (5) Dst, Kp, aa, and other geomagnetic indices; (6) the equatorward edge of the diffuse or discrete aurora; (7) sudden ionospheric disturbances (as inferred from Solar Flare Effects (SFEs); and (8) nitrate deposition in ice. As we move through these parameters, we shift from causes to consequences, and, from a societal point of view, it is those parameters affecting the near-Earth environment that have the most impact. The connection between different parameters also drops significantly as they become further removed from one another. A massive flare on the backside of the Sun, for example, will not likely leave a record in Earth's ice. Additionally, the relationship between parameters need not be linear. Fast CMEs, for example, may not drive an intense geomagnetic storm unless there is a strong negative interplanetary magnetic field (IMF) associated with it.

The term “extreme” is also time scale-dependent. An event might be extreme on the time scale of a year, or decade, but only moderately severe on the scale of 100 years. Here, we invoke a working definition that “extreme” is something we have not experienced within the space era. That is, an event that we have not yet observed during the time our society has become dependent on technology, and which could result in significant adverse consequences affecting a significant fraction of the Earth's population. The Carrington event of 1859 [Carrington, 1859] is perhaps the largest known example of this.

On the morning of 1 September 1859, as Richard C. Carrington was observing sunspots on the solar disk, a particularly large and complex active region destabilized, launching an extremely fast coronal mass ejection toward Earth. A large solar flare ensued, its optical brightness lasting some 5 min and equaling that of the background Sun. Extrapolating knowledge from less severe events observed during the space age, we can infer some of the likely properties of the so-called “Carrington” event. In soft X-rays, the flare has been estimated to be >X10 [Cliver and Svalgaard, 2004]. The ejecta propagated rapidly away from the Sun, generating a fast-mode wave ahead of it, which rapidly steepened into a fast-mode forward shock. The shock, traveling in excess of 2000 km s⁻¹ [Cliver et al., 1990] accelerated suprathermal ions in the ambient solar wind to high energies. As these accelerated particles streamed away from the shock, they excited plasma waves that pitch angle scattered the ions, further accelerating the particles. Some of these energetic ions escaped, traveling ahead of the shock. As they streamed through the heliosphere, they amplified the ambient resonant plasma waves, simultaneously undergoing pitch angle scattering by them. Propagating through an ever-weakening magnetic field, the particles were focused and decelerated. The first particles arrived at the Earth within an hour, although the peak intensity of the particle distribution arrived with the shock, some 17.6 h later [Cliver and Svalgaard, 2004]. The CME and its associated disturbance rammed into the Earth's magnetosphere, generating one of the largest magnetic storms in recorded history [Tsurutani et al., 2003]. Meanwhile, the energetic particles entered the Earth's magnetospheric environment directly through the polar cap region. The most energetic particles penetrated to the stratosphere and produced nitrogen oxides (NO_x) via impacts with molecular nitrogen and oxygen. The lower-energy particles, on the other hand generated NO_x in the mesosphere and thermosphere. NO_x, acting as a catalyst, destroys ozone in the stratosphere, while HO_x rapidly destroys ozone in the mesosphere, and levels decreased substantially immediately following the arrival of the energetic particles. However, more remarkably, over a period of 4 months, vertical winds enhanced by the polar vortex transported the remaining NO_x from the middle and upper atmosphere into the stratosphere, leading to a second significant depletion in the main ozone layer. Ice core measurements indicated that most of the nitrates were deposited within weeks of production, suggesting that gravitational precipitation (i.e., snow) must have played a key role in its downward transport [McCracken et al., 2001]. While this only partially substantiated description of the Carrington event of 1859 sounds plausible, it is not without controversy. The connection between ice core records and SPEs, for example, is currently under debate [Wolff et al., 2008].

Occurring after important advances in observing techniques and theory had occurred, but before the advent of many more, the 1859 storm holds the record as the largest space weather event in over 400 years [McCracken et al., 2001]. Additionally, since the event occurred only ∼150 years ago, it is a constant reminder that a similar event could reoccur any day. But how probable is such an event? With only one such case in the historical records, it cannot be readily predicted from simple “time to event” calculations.

Several recent studies have begun to address the probability of extreme space weather events. Ruzmaikin et al. [2011] studied the power law distribution of coronal mass ejection (CME) speeds, which are an integral component of space weather phenomena. They found that, within the range of 700–2000 km s⁻¹, the speeds were distributed according to a power law with exponent, α = −3.4. They argued that the existence of such a self-similar distribution suggests a single physical process for their generation (at least within this speed range). Additionally, they noted that the time intervals between the events are not independent: Fast CMEs are produced in clusters (a result previously inferred from observations by Feynman [1997]).

Barnard et al. [2011] inferred a rate for major SEP events of P_SEP = 5.2 per century, based on nitrate measurements from ice cores that contained 15 major events between 1700 and 1970 [Shea et al., 2006]. Restricting their interval to the space era, however, which has thus far produced only one major event led to a probability of occurrence, P_SEP = 2.6 per century. McCracken and Beer [2007] also commented on this apparent decline using ground-level event (GLE) data. However, whether this represents an underlying trend, or merely the result of low-number statistics remains to be determined.

A number of techniques can, in principle, be applied for estimating the likelihood of a rare event. These include: event trees, similarity judgments, “time to event,” and extrapolation from more frequent events. Event trees rely on breaking up a chain of smaller events that lead to the catastrophic event being studied, with each smaller event being connected to another either by an “and” or “or” gate [Wang and Roush, 2000]. The probability of the main event is then determined by a combination of the probabilities of the smaller events. The inquiry into the Challenger accident represents perhaps the best known example of this type of analysis [Stamatelatos, 2002]. Unfortunately, for space weather applications we do not know the probabilities of the smaller events or their relationship to one another any better than for the main event. Similarity judgements [e.g., Bordley, 2011] extend the probability of known rare events to new situations, such as estimating the probability of an earthquake on the east coast based on more available data from earthquakes on the west coast, say. Again, however, for space weather situations, there are no obvious known events that can be linked to the one under study. “Time to event” techniques use the time between the occurrence of two or more rare events to estimate the probability of an event happening. Unfortunately, for events that are rarely, if ever observed, estimating this time may not be possible. Even so, such calculations may provide a useful basic check on results obtained from other approaches.

Extrapolation techniques, unlike event trees or similarity judgments, do not rely on knowing any underlying physical processes giving rise to the rare event, or the availability of two or more case studies. Instead, they rely on the assumption that the range over which the events are well observed can be reliably extended to regimes where they are rarely, if ever observed. In this study, we will argue that many space weather parameters exhibit a power law distribution that can be extrapolated to extreme regimes. With this assumption, it is straightforward to estimate probabilities.

Finally, we note that there is a useful distinction worth making between predicting specific rare events and assessing the probability that one will occur within a certain period of time. For the types of space weather phenomena that concern us here, the former task is incredibly more complex and we will focus on probabilistic assessments from hereon.

2. Methodology

Here we outline the basic tools we will employ to compute the probability of occurrence of an extreme space weather event. A set of events, x, is said to follow a power law distribution if the probability of occurrence, p(x), obeys the following relationship:

where the exponent, α, is some fixed value, and C is a constant determined from where the power law intercepts the y axis. Known as Zipf's law or the Pareto distribution, power law distributions, which relate the magnitude of an event to its frequency, appear remarkably frequently throughout all areas of science, including earthquake magnitudes [Christensen et al., 2002] and solar flare size [Lu and Hamilton, 1991]. It is worth emphasizing that power laws fall off much less rapidly than the more often encountered Gaussian distribution. Thus, extreme events following a power law tend to occur far more frequently than we might intuitively expect [Newman, 2005].

Following McMorrow [2009], a more useful quantity for our study, however, is the complementary cumulative distribution function (CCDF), P(x), which is defined as the probability of an event of magnitude equal to or greater than some critical value x_crit:

which, using equation (1) becomes

Thus, the CCDF also follows a power law with a reduced exponent (α − 1).

There are several advantages in using CCDFs over the original power law distributions. First, they avoid the problem noisy tails. Unless p(x) is computed from bins of data that are logarithmically spaced, the ever-smaller number of events landing in the largest event bins drives the errors associated with those estimates up [Newman, 2005]. On the other hand, when the data are summed in a complementary manner as with the CCDF, the noise is minimized. In fact, the data need (and should) not be binned at all; the CCDF at some value x is just the sum of the number of events where x ≥ x_crit. Second, computing the slope of the power law distribution (which is the single most important parameter describing the data) from a least squares fit is prone to systematic biases [Goldstein et al., 2004]. Instead, the maximum likelihood estimate (MLE) is a significantly more accurate method for computing the slope:

where x_i are the measured values of x, N is the number of events in the data set, and x_min is some appropriate minimum value of x, below which the power law relationship breaks down [Newman, 2005].

A third advantage of the CCDF relates to its application. When considering extreme space weather events, we would usually like to know the probability of occurrence of some event of a particular strength or greater. That is, the consequences of an event larger than a certain threshold, such as the Carrington event, are what concern us from a societal view, rather than the probability of an event say between x and x + Δx.

Given the probability of an event as large as or larger than x_crit it is simple to compute the number of events expected to occur during the period covered by the data set:

where N is the total number of events in the data set. Moreover, again, assuming the events occur independently of one another, we can use the Poisson distribution to infer the probability of one or more events greater than x_crit occurring during some time Δt:

where τ is the total time span of the data set. Equations (3)–(6) thus give us a robust method for computing the probability that an event of severity exceeding x_crit will occur some time within the next Δt years, subject to several assumptions.

The assumption that the data follow a power law, both within the regime where measurements are made and into the region of rare and, for the most part, unseen events must be carefully considered. Typically, the power law relationship must break down at both low- and high-frequency events. In the high-frequency portion of the spectrum, the curve usually flattens often due to the fact that smaller events are less easily measured or identified. At the low-frequency portion of the spectrum, several factors could play a role. First, the statistics of small numbers may lead to the curve deviating from what might otherwise be a straight line. Second, and particularly when the power law distribution falls off more rapidly than higher-frequency rates would suggest, the cutoff may be a real physical limitation. For example, there is undoubtedly a limit to the speed of CMEs, related to the maximum kinetic energy that can be derived from the finite magnetic fields within an active region. Even if arguments could be made to circumvent this, as an absolute limit, CMEs certainly cannot travel faster than the speed of light. Third, the distribution may be composed of distinct parts, each governed by different physical processes, such that the curve becomes either steeper or more shallow than the main portion that is measured. The initiation of coronal mass ejections, for example, may occur through flux cancelation, “break out,” the kink instability, flux emergence, or some other mechanism [Forbes, 2000], each of which could produce CMEs with different distributions of speeds. Similarly, energetic particles are thought to be produced either through an impulsive mechanism or a gradual one [e.g., Lin, 2011].

A second implicit assumption in our analysis is that of time stationarity; both during the time of the measured data, and over the time in the future over which we aim to extrapolate. For the Sun, this is clearly an approximation that requires careful assessment. On the scale of a decade or so, the solar activity cycle modulates many solar parameters [e.g., Riley et al., 2000]. The largest 2% of geomagnetic storms (the so-called “super storms”), for example, tend to occur shortly after the maxima of the decadal-scale solar activity cycles, occurring most often near the equinoxes [Bell et al., 1997]. Thus, the data used to construct the power law should have a duration at least as long as this. By extension, we must also recognize that any predictions made over say the next decade would necessarily be solar cycle-averaged predictions, and the actual probability of an extreme space weather event at some point in the cycle may be different.

A final tool that will be useful for our analysis relies on the average time to the next event to compute a probability of occurrence. For Bernoulli distributions, that is, independent events that either happen or do not, with a constant probability of occurrence, it can be shown that the probability of occurrence is given by

where τ is the average time to the event. Thus, an event that occurs once every 100 years would have a probability, P = 1/(1 + 100/10) = 0.09, or 9% of occurring during the next decade.

3. Analysis of the Data Sets

Although there is an almost endless list of candidates, we have chosen four space physics data sets from which to assess the likelihood of an extreme event with respect to those measurements: Peak rates from solar flares; the speeds of CMEs; the strength of geomagnetic storms, as determined from Dst; and nitrate spikes in ice cores, arguably capturing large SPEs. These choices encompass a range of “causes” and “consequences” along the chain from the solar surface to the polar ice. And, while their relationship has not been well established, flares and CMEs are clearly related to the same explosive release of magnetic energy in the corona [Gosling et al., 1993]. Similarly, SPEs likely represent a subset of SEPs. Our final motivation for choosing this subset of parameters was to highlight both the power and the limitations in applying this power law extrapolation technique.

3.1. Solar Flares

X-ray fluxes from solar flares represent a natural starting point for analyzing the probability of occurrence of extreme space weather events since it has been firmly established that they follow a power law distribution in peak photon flux over many orders of magnitude [Lin et al., 1984; Dennis, 1985; Lu and Hamilton, 1991].

For our analysis, we use the hard X-ray measurements from the BATSE instrument on board CGRO (ftp://umbra.nascom.nasa.gov/pub/batse/), which recorded solar flare data continually for almost a full decade. In Figure 1 we show the peak count rates (proportional to flux) of all 7236 BATSE X-ray flares observed from April 1991 through May 2000. We note several points. First, the time series is clearly nonstationary: During the approximately solar cycle duration of the data, the peak rate decreases to a minimum in mid/late 1996 and then rises, in phase with sunspot number. Second, the events are strongly biased toward low peak rates, with only a tiny fraction exceeding 10⁵ counts/s/2000 cm².

Figure 2a shows a histogram of the number of events as a function of peak rate. Several points are worth noting. First, the noise increases as the peak rate increases. These data were binned in intervals separated equally in linear, not logarithmic space, and thus the number of events at high values of peak rate becomes extremely small (→ 0). Although this could have been partially rectified by binning the data in logarithmic space, it serves to highlight one of the advantages of computing the CCDF: The noise does not increase with increasing peak rate since all of the data at a particular peak rate and above are summed to compute that value. Second, although the histogram shows a clear power law distribution over several orders of magnitude, at lower values (<10⁴ counts/s/2000 cm²) the slope becomes shallower. Furthermore, below a peak rate of ∼5 × 10² counts/s/2000 cm², the number of events drops quickly to zero.

In Figure 2b we show the CCDF, i.e., the probability of an event occurring that exceeds some critical peak rate, as a function of peak rate. The noise is significantly reduced. Here, the change in slope is also apparent, as is the rollover at 10² s/2000 cm², which appears here as a flattening in the curve. We note that the slope of the event-frequency plot is ∼−1.8, in agreement with previous studies [Lin et al., 1984; Dennis, 1985], while the slope (computed using MLE) for the CCDF is −0.84. Theoretically, we would expect the latter to be one less than the former, which, given the errors associated with computing the former, is relatively consistent. More importantly, however, for slopes where (α − 1) < 1 (or α < 2), the mean value of the peak rate becomes infinite [Newman, 2005]. Of course the BATSE data contains a finite number of samples, and the mean peak rate must therefore remain finite. In practical terms, this suggests that mean values computed from a number of subsets within the BATSE data set would not converge.

The slope of the CCDF has two major issues for predicting future extreme events. First, there appears to be either two slopes or a general curve from a more shallow to steeper slope over the well-populated portion of the data (from 10³ to 3 × 10⁵ counts/s/2000 cm²). We used the second portion of the distribution to fit a straight line (using both least squares (dashed) and MLE (solid). However, more significantly, the data drop sharply at higher peak rates, suggesting a saturation in flare peak rate at ∼2 × 10⁶ counts/s/2000 cm². While the theoretical implications of this are interesting to consider, for our purposes, this proves to be a significant hinderance for estimating the probability of extreme events. Although it is possible that this is the result of errors stemming from low-number statistics (which we will argue is likely the case for another data set later), here, the deviation from the power law line begins while there are still a reasonably large number of events and the profile is smooth.

The next issue we encounter lies in our choice of a critical peak rate from which to compute a probability of occurrence. For the 1859 event, we have only indirect evidence about the size of the solar flare observed by Lord Carrington. In fact, he is credited with the first observation of a solar flare. Unfortunately, Carrington observed it in white light and described it only in qualitative terms [Carrington, 1859]. And, although it must have been extremely strong to have been visible to the naked eye, Cliver and Svalgaard [2004] estimated that it was only among the top ∼100 flares of the last century and a half, undoubtedly exceeding X10.

Ideally, we would like to use measurements of brightening in visible light, which is what Lord Carrington observed [Carrington, 1859]. However, so-called “white light” flares are too rare to generate accurate distributions. Hard X-ray emissions, on the other hand, have been regularly observed for decades. This emission, though, is produced by nonthermal electrons in the chromosphere, as opposed to the white light emission, which originates near the photosphere. Watanabe et al. [2010] found that the energy in white light emission was equivalent to the energy supplied by all electrons accelerated above 40 keV, suggesting that the electrons are responsible for the white light emission. However, beyond this, translating what Carrington qualitatively described into a robust peak flux rate in X-rays is challenging, if not impossible.

Although we could apply equations (3)–(6) to estimate the probability of a large event occurring, given that (1) we have no reliable “yardstick” for how big the Carrington event would have been in hard X-ray emission and (2) the power law relationship appears to break down at peak rates exceeding 2 × 10⁶ counts/s/2000 cm², any probabilities derived would be dominated by uncertainty.

3.2. Speeds of Coronal Mass Ejections

The strongest geomagnetic storms are produced by fast CMEs, and faster CMEs produce more severe effects [Gosling et al., 1990]. Thus, a next logical stop in the Sun-to-Earth chain is to consider the distribution of CME speeds. For this, we have extracted speeds from the LASCO CME database (http://cdaw.gsfc.nasa.gov/CME_list/), which contains plane-of-sky speeds for 14,735 CMEs observed over the lifetime of the SOHO mission (1996–2010). For this study, we use the quadratic speeds, which are computed by fitting a second-order polynomial fit to the height-time measurements and estimating the speed at the time of the highest measurement possible. In Figure 3 we show CME speeds as a function of time. As with the X-ray flare data, we note the strong clustering at low speeds. Additionally, there is a trend for the number and speed of the fast CMEs to increase from 1996 (solar minimum) to 2002–2004 (2 years after solar maximum) returning to a relative minimum again in 2008–2009 [Riley et al., 2006].

Cliver et al. [1990] estimated that the 1859 event took 17.6 h to travel the distance from the Sun to Earth, suggesting an average speed of ∼1.5 × 10⁸ km/(17.6 h × 3600 s/h) = 2400 km s⁻¹. In comparison, from the LASCO database, two of the fastest CMEs of solar cycle 23 (on 10 April 2001 and 24 September 2001) had travel times of ∼32.5 h. Since both of these events had observed plane-of-sky speeds of ∼2400 km s⁻¹, to a first approximation, we might estimate an initial speed of ∼2400 × 32.5/17.6 = 4500 km s⁻¹ for the 1859 CME. Since an initially faster CME will undoubtedly display a stronger declaration profile than this linear extrapolation presupposes, we could reasonably round our guess for the Carrington CME's initial speed to 5,000 km s⁻¹.

In Figures 4a and 4b we show a histogram and CCDF of CME speeds. We note the relatively good power law distribution in speeds in Figure 4a as well as the rollover, marking the median speed of events, at ∼350 km s⁻¹. In agreement with the study by Ruzmaikin et al. [2011], we find power law spectra both in the histogram and CCDF over a relatively broad range between 700 km s⁻¹ and 2,000 km s⁻¹. Using this slope (−3.2), for a critical CME speed, v_CME ∼ 5,000 km s⁻¹, this would suggest a probability of observing such an event, or greater, over the next decade as 85%. However, the fastest CME observed by LASCO in 15 years of operation was only 3500 km s⁻¹, and thus, we believe this estimate is not credible.

We suggest that the reason for this artificially high probability is that above 2000 km s⁻¹ there appears to be a well-defined “knee” in the distribution. Therefore, to address this, in Figure 5, we have replotted the CCDF for the highest speeds and computed a MLE fit to only the distribution above 2000 km s⁻¹. Thus, using the more restricted fit to the tail of the CME speed data, the same extreme event, (i.e., v_CME > 5,000 km s⁻¹) would only have a probability of 12%, based on a revised slope of −6.1. (For reference, we note that the least squares fit to these data produced a slope of −6.5, the most significant difference of all four data sets, leading to a probability estimate of 8.5%).

3.3. Geomagnetic Storms

In general terms, a geomagnetic storm is a disturbance in the Earth's magnetosphere driven by changes in the solar wind. Both high-speed streams as well as CMEs can drive such storms, with the latter typically leading to more significant effects. Different components in the solar wind produce different consequences in the magnetosphere. Fast-mode shocks and their following sheath, for example, compress the magnetosphere, while pronounced intervals of southward interplanetary magnetic field can transfer energy from the solar wind to the magnetosphere [Kivelson and Russell, 1995].

Although there are a number of indices that capture various aspects of geomagnetic activity, large negative excursions of the “disturbance–storm time” index, or Dst, perhaps best describes the main phases of a magnetic storm [Gonzalez et al., 1994]. Moreover, it is arguably the best “societal impact” parameter. Although anecdotal, we remark that while the 13 March 1989 event, with a peak Dst < −600 nT, caused the collapse of the Hydro-Québec power grid, and a resulting loss of electrical power to six million people [Bolduc, 2002], the so-called “Bastille Day” event of 14 July 2000 which was associated with a peak Dst of −300 nT, caused no power failures or other significant terrestrial effects. Tsubouchi and Omura [2007] used Dst measurements between 1957 and 2001 to estimate the probability of occurrence of intense geomagnetic storms, finding that a storm on the scale of the 1989 event was likely to occur every 60 years or so.

A “storm” can be arbitrarily defined when the main phase falls below some value, typically −50 nT. Additionally, we can classify storms as moderate (−50 nT > Dst > −100 nT), intense (−100 nT > Dst > −250 nT), and severe (−250 nT > Dst > −600 nT). The Carrington event of 1859 was initially estimated to have a peak negative Dst of −1760 nT [Lakhina et al., 2005], although this was later revised, and reduced by a factor of two to −850 nT [Siscoe et al., 2006]. Since the beginning of the space age (1958), only one storm has exceeded −600 nT, the 1989 storm for which Dst ∼ −640 nT [Gonzalez et al., 1994].

In Figure 6 we show a 46 year time series of Dst. Several points warrant comment. First, Dst is strongly asymmetric: storms reveal themselves as sharp negative excursions lasting a few days at most, and hence, on this scale, a single vertical line. This allows us to relax the need to track the sign of Dst since whenever |Dst| exceeds, say, 100 nT, it must be a negative quantity. Second, Dst's envelope is clearly modulated on a time scale of ∼11 years, with minima (maxima) in the amplitude coinciding with solar minima (maxima).

To generate an event-based data set, we arbitrarily define a “significant” magnetic storm to be one for which |Dst| exceeds 100 nT. In principle, we could define an “event” as the hourly value of Dst and compute our probabilities based on that. However, this might impact the results in two negative ways. First, we are not concerned with the slope of the CCDF for small values of Dst. Therefore, this portion of the parameter space is not relevant for computing the slope. Second, we are interested in predicting the probability of events, where one event is when Dst exceeds a threshold for some period of hours or days. Thus, we would rather identify a contiguous range of data that all exceeds some criteria as a single event, rather than a set of events. In Figure 7, we show the occurrence of these significant storms as a function of time. Again, their density distribution changes in phase with the solar cycle. We note that only one event exceeds 600 nT, and, moreover, only two events exceed 400 nT. Another potentially important trend is the tendency for the strongest storms to become stronger over the last four solar cycle maxima. Whereas the top 5 storms around 1970 lay between 200 and 300 nT, by 2005, the five most intense storms lay between 350 and 450 nT. If such a trend is real it suggests a basic nonstationarity in the data on the same scale as the duration of the data set, and that predictions of future events may underestimate their true probability.

In Figure 8a, we show a histogram of events as a function of the severity of the storm. The data appear to follow a power law distribution, as indicated by the least squares fit to the data. In Figure 8b the power law relationship of the CCDF is considerably better: Only the last 3 points (which are made up of only 1, 2, and 3 events, respectively) deviate. The slope of the MLE fit is −3.2. Again, we note the clear advantage in using the CCDF as computed here, rather than a histogram approach. The histogram was computed using 100 bins equally spaced in linear Dst space, while the CCDF on the left was created by cumulatively summing points from the right. Therefore, the CCDF contains all original 746 geomagnetic storms that were identified in the data set, whereas the histogram contains only 100 points.

Using equation (3), we compute the probability of an event of magnitude equal to, or greater than a threshold value of Dst_crit = −850 nT to be 0.001. And, using equations (5) and (6) we estimate the probability of such an event occurring over the next decade to be 12%. If we require an event to exceed a threshold of −1700 nT, the probability of it occurring over the next decade reduces to 1.5%.

3.4. Extreme Space Weather Events From Nitrate Records

Finally, in the chain of space weather parameters from the Sun to the Earth, we arrive at space weather records potentially contained within ice cores. The value of these data lies in their long time span, going back more than 400 years; however, they are not without caveats. First, while the nitrate spikes are generally believed by space physicists to be a record of large, historical space weather events [McCracken et al., 2001], ice core chemists are skeptical [Wolff et al., 2008]. They posit that no viable mechanism exists by which Solar Proton Events could be imprinted within the ice, suggesting instead that high concentrations of sea salt provide a simpler and more consistent explanation for the deposition of aerosol nitrates. Second, there are only 70 events spanning the 450 years for which we have data. The largest event in the data set, with a fluence of 18.8 × 10⁹ cm⁻¹, occurred in 1859. That is, the largest event in the last 400 years was the Carrington event. More importantly, however, with such a limited number of events, the statistics of the fit and the resulting probability estimates will be more prone to error.

Figure 9 shows all >30 MeV proton events with fluences exceeding 10⁹ pr cm⁻² as a function of time between 1562 and 1944. Although it is not possible to show rigorously, because of the limited sample size, there is no obvious trend in the distribution of event sizes or temporal clustering to suggest that the time series is obviously nontime stationary. Of note is that the Carrington event is substantially larger than the other events in the data set, with the second largest event producing a fluence of only 59% of the value of the 1859 event.

Figure 10a shows a histogram of events versus event size and Figure 10b shows the CCDF as a function of event size. Comparison of the two distributions emphasizes the strength of the CCDF approach. Whereas 70 events contributed to the plot on the right, the plot on the left was made up of only 9 bins, a number large enough to see a trend in the data, but small enough that a sufficient number of points would fall into most of the bins. The MLE fit to the line in Figure 10b gives a slope of −2.0. Using equations (3)–(6) we estimate the probability of an event of, or exceeding 18.8 × 10⁹ cm⁻¹ occurring over the next decade to be 3.0%.

4. Summary and Discussion

In this study, we have applied a power law probabilistic analysis to assess the likelihood of a space weather event on the scale of, or larger than the Carrington event of 1859. Based on Dst variations over the last 45 years or so, we inferred a probability of ∼12% that an event with |Dst| > 850 nT would occur over the next decade. By requiring a more significant threshold of 1700 nT, the estimate reduced to 1.1%. Similar analysis of CME speeds also yielded a probability of 12%. For nitrate records, for which the Carrington event is believed to be contained within, a probability estimate of 3.0% was found. We also investigated hard X-ray flare data, but, because of several significant limitations, we were unable to obtain a reliable estimate.

Although it is tempting to ascribe extra significance to the matching probabilities we obtained from Dst and CME speeds (12%), we stress that this may be more of a coincidence than any underlying truth. It is possible, for example, that the CME responsible for the Carrington event was only traveling at 4,500 km s⁻¹ and that Dst ∼ −1000 nT, in which case the probabilities would have become 21% and 7.3%, respectively. Nevertheless, based on our initial estimates for Dst and CME speed, a probability of ∼12% remains our best estimate.

Our analysis has relied on several important assumptions, the most significant of which is that the data are time stationary. On the shortest time scales of a decade or so, the solar cycle drives important variations in most, if not all of the quantities we studied here. Over longer time scales, there is ample evidence for nonstationarity. Solanki et al. [2004], for example, argued that compared to activity over the last 11,000 years, the last 70 years have been a time of exceptional activity. Steinhilber et al. [2008] also found that solar activity now is stronger than 85% of the time over the last 9,300 years. Nonstationarity can also be seen ‘directly’ in ancient proxy sunspot number catalogues [Wittmann and Xu, 1987]. If the time series were not stationary in the past, they are no more likely to be stationary in the future. With a consensus beginning to emerge between various reconstructions of the heliospheric field on the time scale of several centuries [Svalgaard and Cliver, 2010; Lockwood and Owens, 2011], it may be possible to predict the future large-scale behavior of solar activity. Abreu et al. [2008], for example, estimated that the current grand maximum is only likely to persist for another 15–36 years. Similarly, Lockwood et al. [2009] argued that solar activity rose during the first half of the 20th century, peaking in years 1955 and 1986, and subsequently declined, with the grand maximum ending within the next 20 years or so. Lockwood et al. [2009] further estimated that there is a 10% probability that activity on the Sun will decrease producing grand minimum conditions during the next 40 years. Obviously, if such conditions do ensue, probabilistic forecasts based more active solar conditions may be less accurate. However, not necessarily for society's benefit: During periods of very low activity, for example, radiation from galactic cosmic rays (GCRs) will be higher, posing larger risks for passengers and airline crew, as well as avionics. And, while SEP events themselves may decrease, the consequences from the ones that are produced may be more severe [e.g., Barnard et al., 2011].

A related but distinct effect is clustering. CMEs have been shown to occur in groups [Feynman, 1997; Ruzmaikin et al., 2011], which, to some degree, violates the assumption of time stationarity. Moreover, from a physical perspective, the clustering of events may lead to interactions between them as one CME overtakes another ahead of it. This may destroy or dampen a previously extreme event, or, conversely, lead to a nonlinear amplification, promoting a large event to an extreme one.

A second important assumption for our analysis is that the power law relationship extends beyond the measurements, that is, that events more extreme than have been observed continue to fall along a straight line in log-log space. We argued that the flare data (Figure 2) deviates significantly from the power law relationship that holds for lower peak rates, and, thus, cannot be reliably extrapolated beyond observed events. In contrast, the remaining quantities could, arguably, be considered to lie on power law lines up to the most severe observed case. Moreover, the solar proton events, which are the only measurements that captured an extreme event (the Carrington event), displayed power law characteristics at all fluences.

Each of the parameters we chose to analyze in this study displayed unique strengths and weaknesses for the purpose of computing a probability estimate of extreme space weather events. The nitrate record, which is the only one that actually purports to have captured the Carrington event is limited to a small number of events, making statistical inferences the most dubious. However, issues have been raised about whether the ice core record is actually measuring events from space [Wolff et al., 2008]. In contrast, the hard X-ray data from CGRO's BATSE instrument are abundant, but the resulting power law distribution is not of a single slope, making extrapolation dubious. CME speeds, as determined from the SOHO/LASCO database, appear to show two distinct power law distributions, joined at a “knee.” While this may be hinting at the possibility of two classes of CMEs (“fast” and “very fast”), it is also possible that the fastest events are merely an asymptotic falloff. We currently lack the understanding of CME initiation and acceleration to be able to place robust limits on initial CME speeds. Finally, analysis of Dst perhaps yields the best estimate for the prediction of the likelihood of an extreme event. It is global in scope, contains a large number of events, spans four solar cycles, and appears to follow a single power law distribution. However, even here, we must reiterate one of the primary assumptions of this analysis; that the power law can be reliably extrapolated beyond known events. However, here, we also have indirect estimates of Dst during the Carrington event that provide at least a qualitative check on our results.

It is also not surprising that the probabilities we have estimated for another Carrington event occurring within the next decade based on Dst and nitrate records do not match: In addition to the aforementioned controversy about whether the nitrates records even capture space weather events, they are predictions of different consequences.

Our choice of data sets and the probabilities determined should not be construed as definitive, in any way. The parameters were chosen to illustrate a range of distributions, features, and limitations. Instead, we could have analyzed soft X-ray data from GOES [e.g., Feldman et al., 1997], the Solar Proton Event (SPE) list maintained by NOAA (http://www.swpc.noaa.gov/ftpdir/indices/SPE.txt), or the equatorward edge of the auroral oval [e.g., Newell and Meng, 1988]. For each data set, we would have identified similar, but unique artifacts and derived similar, but distinct probabilities.

The choice of the best parameter from which to compute these probabilistic predictions depends on several factors. First, they should be closely related to the effect one is trying to predict. SPEs, for example, would be an appropriate parameter for estimating radiation doses of astronauts, while localized measurements of the horizontal component of the magnetic field might be more suited for ground induced currents (GICs) and their effects on the power grid. Second, the parameter should be distributed according to a power law. Although solar flares are perhaps the best known example of a power law distribution in space physics, we showed that at least the peak rate from BATSE measurements is not ideal. Third, one should not be able to present a credible argument that the power law relationship breaks down before the critical value to be predicted. It may be that the best parameter is not something that is usually studied. For example, it is well known that geomagnetic activity is sensitive to the speed of the CME, but also the magnitude and sense of B_z, i.e., the dawn-dusk electric field across the magnetosphere. Therefore, an alternative, and arguably better parameter that could be studied and predicted would be ‘events’ where the average solar wind electric field drops below some threshold for some minimum period of time, provided that these ‘electric field events’ follow a power law distribution.

Risk analysis involves balancing the probability of a particular event occurring with the consequences of that event. While we have focused on the first part of this relationship, others have speculated on the consequences that an event at least as large as the Carrington event could have on society [Schieb, 2011]. Unfortunately, the predicted effects are at least as poorly known as the probability of occurrence, thus limiting the ability of decision makers to use these results to guide policy. Nevertheless, it may be worth putting the probabilities computed here into perspective. Consider, for example, the probability that an asteroid (of sufficient size that it would cause substantial damage to a localized region of the globe) will hit the Earth in the next decade. It has been estimated that objects as large as, or larger than 50 m hit the Earth roughly once every thousand years [Chapman, 1994]. The Tunguska event of 1908, for example, was about this size, and would have been capable of destroying a large metropolitan area [Longo, 2007]. Of course, the average effect of such an event is considerably lower since there is a much larger probability of it hitting an ocean where the death toll could easily be zero. On the other hand, a space weather event will always be global in nature (although its effects would be amplified in certain areas). While global satellite monitoring can now rule out potential impacts in the near future, using equation (7), we can estimate the probability of such an impact over the next decade to be ∼1%, which is an order of magnitude smaller than our estimate for a Carrington event based on Dst and CME speed considerations. The two scenarios become comparable from a risk analysis perspective if the consequences of a Carrington event are even one tenth those of a meteor impact. If financial resources are limited, such knowledge could be used to allocate funds more effectively. Moreover, if better mitigation strategies exist for one scenario than another, they could be more aggressively pursued. For example, if major, prolonged power outages, which could, in turn, result in food and water shortages, potentially over many months, are realistic consequences of a Carrington event [Schieb, 2011], effort should be expended into building infrastructure that can withstand such an event. In contrast, there appears little that can be done, at least in the near term, to mitigate a meteor impact, beyond effective evacuation measures.

Computing probabilities for extreme space weather events may ultimately be useful in shaping policy decisions. One of the problems with extreme events is that prior to their occurrence, their perceived risk is effectively zero, yet following it, the risk rises to nearly 100%. Thus, it is important to develop objective estimates for these probabilities, which can serve as a starting point for risk analysis, even if the errors associated with them are larger than we might like. Thus, if the results presented here are further substantiated, are there measures that could, or should be enacted to minimize damage from an extreme solar event?

In closing, we reiterate that our primary aim in this study was to introduce a technique for estimating the probability of occurrence of extreme space weather events. Additionally, our analysis has shown that a relatively rich subset of space physics data can be approximated by power law distributions. Our results allowed us to answer a basic question, at least in an approximate way: How likely are such events? The answer of course depends on what you mean by “event” and how severe you define “extreme” to be. Nevertheless, our results overall suggest that the likelihood of another Carrington event occurring within the next decade is ∼12%.