TESTING AUTOMATED SOLAR FLARE FORECASTING WITH 13 YEARS OF MICHELSON DOPPLER IMAGER MAGNETOGRAMS

Awareness of space weather and its implications for sensitive technologies has been steadily increasing. As society becomes more technologically dependent on complex global systems, the potential risk posed by the geomagnetic response to solar variability increases. Adverse effects include large-scale power distribution failure, Global Positioning Satellite tracking errors, high-frequency radio communication interference, and low-orbit satellite drag (NRC 2008). These issues originate in the interaction between the interplanetary magnetic field in the solar wind and Earth's magnetosphere, as well as powerful radiation from solar flares. The spawning point for these phenomena is almost always solar active regions (ARs).

Schrijver (2009) described a process of solar flare initiation based on previously published models. When field lines snap to a lower energy configuration via reconnection, the excess energy released produces coronal mass ejections (CMEs) and/or solar flares (e.g., Silva et al. 1996; Metcalf et al. 1995). Nonpotential energy is stored in the magnetic fields (e.g., Tanaka & Nakagawa 1973; Krall et al. 1982), parameterization of which is an active area of research (e.g., Falconer et al. 2002, 2003, 2008; Leka & Barnes 2003a, 2003b; Barnes & Leka 2006). Schrijver (2009) also calls for more statistically significant studies based on consistent and long-duration observations of ARs, such as the present analysis.

This investigation explores the magnetic field parameterizations proposed in previous studies. To date, there are no consistent, long-duration vector-field observations available. Thus, the parameterizations we chose to study were restricted to quantities derivable from longitudinal magnetograms. In addition, we wanted parameters robust enough to be calculated with little or no user interaction in order to (1) facilitate the application to a large data sample and (2) remove user-generated bias. After obtaining preliminary results, parameters that lacked an indication for further analysis were no longer pursued (see Section 2.2.4). The most promising parameters that satisfied the above criteria were the total unsigned flux, bipolar region separation, and the gradient-weighted inversion-line length (GWILL).

These parameters confirm the findings of other investigations over many years. Guo et al. (2006) tested a possible measure of nonpotentiality deemed the effective distance, the distance between the flux-weighted centers of the bipolar region constituting the AR, and found that it correlates well with flaring intensity. Falconer et al. (2002) tested the correlation between AR CME productivity and several possible measures of nonpotentiality, such as the length of the strong-shear, strong-field primary inversion line. Falconer et al. (2003) generalized the measure for line-of-sight magnetograms. Schrijver (2007) defined a similar measure called R that parameterized the unsigned flux near high-gradient, strong-field polarity-separation lines. Leka & Barnes (2003b) showed that the predictive capability of any of these measures based on longitudinal magnetograms is not particularly strong, even when multiple measures are combined to optimize their predictive ability.

Most previous studies have either been restricted to a relatively small sample size or neglected the time evolution of the magnetic field. Leka & Barnes (2007) showed that removing these two limitations is critical to improving the ability to understand and predict solar flares. This study expands the earlier work to a time-sensitive and statistically significant sample using the 13 year record of Michelson Doppler Imager (MDI) synoptic line-of-sight magnetograms.

2.1. Initial Data Collection and Extraction

This study uses the recalibrated level 1.8, 96 minute cadence line-of-sight full-disk magnetograms from the MDI instrument (Scherrer et al. 1995) on board the Solar and Heliospheric Observatory (SOHO). These products are available from 1996 to the present with relatively few gaps. The present investigation analyzes the 1075 ARs visible in the record that span most of the interval from 1996 April 15 to 2008 December 31.

Very few selection criteria were applied to this initial data set. All of the numbered National Oceanic and Atmospheric Administration (NOAA) ARs within this time period were considered. ARs needed to be observed in at least three magnetograms within 30° of the central meridian to be included. This reduces the influence of projection effects that distort the MDI measurements far from the disk center. Erroneous and incomplete images were removed.

The "FindAR" program described here uses the Solar SoftWare IDL library (Freeland & Handy 1998) and is available online.¹ The only required user input for FindAR are NOAA AR numbers, which can be input individually or by range. FindAR can also accept a date range and compile the relevant NOAA AR numbers automatically. The program accesses the online NOAA daily solar regions database, which provides the location, type, and approximate size for each region. The catalog information is used to identify and extract a rectangular sub-image from each of the full-disk MDI magnetograms that includes the AR. A buffer area was added around each rectangular AR to ensure that all elements of the AR are included. The buffer size scales with the size of the AR and is adjusted to maintain a constant re-mapped image size for a given AR during the disk passage. The MDI group's "project" module was used to extract sub-images and re-map the image using a Postel equidistant azimuthal projection. This map projection yields consistent estimates of distance on the solar surface, which was necessary for the reliable calculation of the magnetic measures. Figure 1 shows an example of a series of sub-images as the AR crosses the disk. Calculations were performed on the remapped images to characterize the evolving AR parameters. These calculations were performed using our "Solignis" program, which includes the code for FindAR and requires no additional user input. Solignis generates all the necessary data to perform a statistical analysis of all the parameters included in the code. This code can be found at the same URL as FindAR and can be used as the framework for other researchers to perform similar long-duration, consistent database studies with the caveat that they must have access to the MDI full-disk magnetograms.

Zoom In Zoom Out Reset image size

2.2. Calculated Parameters

This section details the four parameters chosen as proxy measurements of the nonpotentiality and/or magnetic complexity of the ARs. The selection of parameters was based on their predictive success in previous studies.

2.2.1. Total Unsigned Magnetic Flux

The unsigned total magnetic flux is a straightforward sum of the absolute value of the longitudinal magnetic field measurements. As mentioned in Section 2.1, large ARs have a large buffer area. This effect causes a modest skew in the flux estimate that slightly amplifies the total flux for regions that are already larger (Bokenkamp 2007). However, the field values in the buffer region are typically weak and contribute little to the total flux. Measurements of the strong magnetic field in the umbrae of sunspots (>3 kG) are saturated in the MDI data (Liu et al. 2007) and this tends to bias the total flux value downward. These competing effects are nontrivial to correct and may slightly influence the value of this variable. We apply a minimum threshold of 100 G in our calculation to reduce the effect of noise and to reduce the contribution of the quiet Sun outside the ARs.

Total flux is a commonly used parameter (e.g., Falconer et al. 2002; Leka & Barnes 2003b) that provides a rough measure of the size and strength of an AR. Leka & Barnes (2003b) found that contrary to the conventional wisdom, total flux alone is not a good predictor of flares. We chose to include this parameter to provide additional weight for contradicting the conventional wisdom.

2.2.2. Primary Inversion Line Length

Polarity inversion lines separate patches of positive and negative flux. These occur all over the Sun and are effectively continuous within ARs (see dashed lines in Figures 1 and 2). The primary polarity inversion line (PIL) separates the major polarity regions of an AR. This line is relatively easy to trace by eye but presents a challenge for computer software to quantify. Bokenkamp (2007) developed an IDL program based on a routine written by Falconer et al. (2003) to effectively define the PIL using a three-iteration process. Figure 2 illustrates this process being performed on NOAA AR 10486, which was a delta-class sunspot group.

Zoom In Zoom Out Reset image size

The image is first strongly smoothed and the zero Gauss contours identified. Smoothing effectively removes the small-scale structures that generate extraneous inversion lines away from the PIL. A vector field is then calculated using an alpha = 0 linear force-free field model (Alissandrakis 1981; Falconer et al. 2002). The first stage ends by identifying continuous zero Gauss contours with a horizontal field gradient on both sides above a specified threshold, and with a strong model vector magnetic field strength. This process is repeated for a less smoothed image. The continuous segments of this stage are superimposed on the result of the previous stage (Figure 2, left) and used as the output contour wherever the two iterations agree within a specified number of pixels. Finally, this method is applied to compare the results from a completely unsmoothed image and the results from the previous stage. This procedure emphasizes the connectivity of the line segments, inversion line stability over time, and reduces noisy extensions of the PIL. Inspection of the segment shows that the algorithm produces a fairly accurate representation of the PIL (see Figure 2, right). It is difficult to define a singular PIL in ARs with exceedingly complex geometries, such as delta sunspots (e.g., NOAA AR 10486). In these cases, our algorithm identifies multiple PILs that are a fair match to those that have been subjectively determined.

Conceptually, a longer PIL indicates a more complex magnetic field structure. Furthermore, a large gradient near the inversion line indicates an appreciable difference in the magnetic field over a relatively small distance. This is indicative of shearing or twisting of the magnetic fields. Therefore, we calculate the length, which is a simple sum of the line segments, and the field gradient across the primary polarity inversion line.

2.2.3. Effective Separation

The effective separation was first investigated by Chumak & Chumak (1987). Chumak et al. (2004) and Guo et al. (2006) found the effective distance between the two bipolar regions of an AR to be a useful parameter for flare prediction. First, the two polarities are weighted by the amount of flux they contain and the centers are found. The distance between the flux-weighted centers is indicative of how separated the bipoles are (Guo et al. 2006). A small effective separation is indicative of a compact AR. Such regions are likely to contain a lot of intermingled flux, which may lead to large gradients.

ARs consisting of more complex geometries than a bipolar distribution may produce separations that are difficult to interpret with this type of measure. No modification is applied to the method of calculation for such ARs, which should be considered when weighing the results yielded from this measure. Georgoulis & Rust (2007) proposed a similar measure, B_eff, which effectively describes the amount of magnetic connection within an AR by comparing the net flux to the distance between magnetic footpoints, as determined by the magnetic charge topology model (Barnes et al. 2005). They found a greater than 95% probability of M- or X-class flaring for ARs with critically high values of B_eff.

2.2.4. Gradient-weighted Inversion-line Length

GWILL closely corresponds to the parameter L_SG proposed by Falconer et al. (2003). This parameter is calculated by applying the equation

$$\begin{equation} {\rm GWILL} = \displaystyle \sum \nabla B_{\perp }, \end{equation} \tag{ 1 }$$

where the sum runs along the entire length of the primary polarity inversion line (PIL), and ∇B_⊥ is calculated across the PIL. This measure tends to emphasize regions of the AR that are strongly sheared and magnetically complex.

Schrijver (2007) created a similar parameter, R, that was a measure of the unsigned flux near such high-gradient, strong-field polarity-separation lines. He found a strong correlation between an AR's propensity to flare and the value of log R. Falconer et al. (2003) proposed using the length of the segment of the PIL with a longitudinal field gradient above 50 G per megameter, L_SG, as a measure correlated with CME initiation. L_SG is a proxy for the strong shear length, L_SS, which can only be derived from vector magnetograms. L_SG uses model transverse fields (as GWILL does) in place of the observed vector field used for L_SS. Whereas Falconer et al. (2003) cut off the length of the PIL with the gradient below a threshold, this study simply weights all segments by the gradient. The resultant quantity should be a better predictor of flaring than the standard PIL because GWILL benefits from the best of both the PIL and the gradient across it. GWILL is larger for ARs that have long PILs, which are indicative of a complex field structure, and is also larger when there are high gradients, which is indicative of shearing. Another primary difference from Falconer et al. (2003) is that this study is searching for a correlation with flares instead of CMEs, and thus a direct comparison is difficult.

2.3. Program Outputs

In order to determine a relationship between these indices and the ARs energetic event productivity, we use the online GOES solar X-ray flare catalog available from the National Geophysical Data Center. Each flare has been associated with an AR through their observed locations in Hα images. Flares not associated with any specific AR were not included in our analysis. For each AR, a time series was created for each calculated variable. Figure 3 shows an example with the temporal location of any M- or X-class flares indicated.

Zoom In Zoom Out Reset image size

For each AR, the parameters for every magnetogram included in this study were placed in a single table that is available online.² There are 71,324 magnetograms listed in our Great Table. For each NOAA AR the table lists the set of relevant magnetograms, including the NOAA AR number, date, time, MDI FITS file identifier, Carrington location, pixel coordinates of the sub-image, and the five calculated parameters. Another table listing the date, start time, end time, time of maximum flux, solar location, GOES classification, and the associated NOAA AR number of each of the 12,576 flares is also available.²

The parameters calculated for each of the 1075 ARs are used to search for a relationship to solar flares. Superposed epoch (SPE) plots provide a sensitive tool to discover the average response of a complex, variable system to a certain type of event at a particular set of times (e.g., Chree 1913; Singh & Badruddin 2006). Essentially, the SPE plots produced in this study are the combination of multiple individual time series with the temporal location of a selected set of flares aligned. One can align to any set of key times. The analysis code accepts a list of key times, or a key list, and gathers the available data from the Great Table. For any specific nonpotentiality measure, the time series are aligned relative to the key times to compute the average response to the event (see Figures 4–6). Since MDI takes magnetograms at a 96 minute cadence, our data series exist as a succession of points 96 minutes apart. Separate key lists were constructed for the X-class, M-class, and B–C-class (quiet) flares. The observation of any M- or X-class flares in an AR was cause to exclude that AR from being included in the quiet key list. This reduces the size of the quiet key list, but avoids the problem that larger flares tend to overshadow the effects of the smaller ones.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition, a control set was created that includes the ARs associated with no flaring at all during the ARs passage within 30° of the central meridian (Dead Silent regions). Since an AR could flare while on the far side of the Sun, we chose to enforce the 30° limitation on the Dead Silent definition to keep the analysis consistent with the 30° limitation placed on the selection of MDI data (see Section 2.1). For Dead Silent regions, a random instant in the time series was chosen to act as the key time to eliminate any location-based bias. However, we required that at least three data points be available before and after this random instant (e.g., if the random instant chosen was the first observation within the 30° limit, no prior data would be available and a new random instant must be selected).

If a flaring AR contained more than one flare, we included it multiple times and aligned the time series so that each flare would be at t = 0. In the forecasting community, it is a rule of thumb that an AR that has flared recently is likely to flare again—this is known as flare persistence (Sawyer et al. 1986). At times of particularly frequent flaring, the magnetic fields may not have time to fully relax after a flare before another occurs. This means that the signal for any such flare may be obscured. The analysis considered every flare for an AR instead of choosing a single flare with a predetermined condition. The goal was to improve statistics and detect any weak systematic signal with SPE analysis. The total number of flares in each key list was 98, 1139, and 4932 for the X, M, and quiet key lists, respectively. The total number of ARs lacking any flares was 410, which constituted the Dead Silent key list.

The first four statistical moments of each time in the SPE were calculated. Our primary focus was on the mean and standard deviation. Our goal was to see what tends to happen at each point in time relative to the flare time; thus, we are primarily concerned with the fractional change. Therefore, data shown in Figures 4–6 have been normalized by dividing the time series for each key time by the median value computed for that interval before combining them. We chose to normalize by the median, as opposed to the mean, in order to generate a normalization less sensitive to outliers. We included the Dead Silent regions as a control, as we did not expect that a flaring AR would behave the same as a non-flaring AR. The results confirm this expectation and suggest that the signal is truly a response to the flares.

Figure 4 shows a clear trend in GWILL—a gradual apex centered approximately on the flare time (t = 0) whose magnitude is greater for X-class flares. The control test (lower-right panel) shows that GWILL for ARs lacking any flaring tends to decrease during the disk passage. This decrease may be partially due to a selection bias—regions that are growing are increasing in complexity and are therefore more likely to flare, thus excluding them from the Dead Silent key list. This suggests that on average, an AR's PIL and the gradient across it will decrease when an AR is not producing flares. In Figure 4, note the strong disparity between ARs with flaring and those without flaring. The size of the peak is appreciably larger for X-class flares. Table 1 summarizes the percent changes from −40 and −20 hr to the key time for each quantity.

The effective separation (Figure 5) shows a similar flare dependency to GWILL; however, the general trend is now inverted. The separation definitely decreases prior to flaring and increases very slightly after. The control test (lower right panel) confirms that the separation between bipolar regions remains nearly constant in ARs that do not produce flares.

Finally, the total unsigned flux (Figure 6) demonstrates very similar behavior to GWILL. The primary difference is that the control test (lower right panel) produced a nearly constant value of the total flux. Also, the percent changes for total unsigned flux were markedly smaller than GWILL. Again, X-class flares generate a stronger response than weaker flares. We do not include separate plots for the PIL and the gradient across it because GWILL includes the most responsive parts of both.

The standard deviation for the Dead Silent key list for these three measures is worth noting. In GWILL (Figure 4, lower right panel), the standard deviation is relatively large indicating a lack of signal, as expected. However, standard deviation is not as large in effective separation and total flux (Figures 5 and 6, lower right panel). This difference is likely due to GWILL values that are closer to its fundamental noise level. When values of the parameter are sufficiently low, the standard deviation becomes noticeably larger.

Close inspection of Table 1 will reveal that the percent changes found for the M and quiet key lists are very similar. The quiet key list often has values that are slightly larger than the M key list. This may be due to the greater number of C and B flares compared to M flares. This result suggests that given greater numbers, SPE analysis is sensitive enough to detect a signal, even when the presumed cause of the signal is significantly weaker.

3.1. Forecasting and Prediction

Superposed epoch analysis has the ability to pick out weak systematic responses. While many of the commonly studied very large ARs have been found to have a degree of high-energy-event predictability, it is important to study what happens to a typical active region in the time surrounding such events. In this study, a change in three primary parameters was discovered to be well correlated with the presence of flares. X-class flares are correlated with a stronger response in these parameters. However, the changes are too small to be predictive in an individual case-by-case sense. This is consistent with the findings of Leka & Barnes (2003b), who are skeptical of the possibility of using only line-of-sight magnetic data to predict CMEs or flares. Barnes & Leka (2006) found that the magnetic field at the photosphere is only moderately related to the flare productivity of the region. Barnes & Leka (2008) show that the best measures of magnetic nonpotentiality are as good at forecasting flares as independent methods, and suggest that forecasting may be improved by combining these methods. We concur with these conclusions and urge the forecasting community to carefully consider any prediction results based on longitudinal magnetograms.

In the course of this study, algorithms were developed for the automatic collection and extraction of AR patches from MDI magnetograms, as well as the automatic calculation of variables on all these data. The programs and our results from them have been made available online. However, access to the MDI database is required and editing of the IDL code may be necessary. Development of these routines is underway for use on the Solar Dynamics Observatory's Helioseismic Magnetic Imager (SDO/HMI) data.

We have also begun preliminary research into the change in the total flux over an AR's disk passage. Many researchers have found that rapid changes in flux are well correlated with an ARs propensity to flare (e.g., Wang et al. 1994; Choudhary et al. 1998). We hope to conduct a study similar to the present investigation to analyze this quantity statistically. Petrie & Sudol (2009) have shown that local small-scale changes within the AR are associated with flaring. This suggests that some information may be lost when parameterizing an AR image with a single number. Future statistical studies should consider including this type of spatial analysis. Many of the most promising flare-predictive magnetic parameters require the use of vector data, which are not currently available in a long duration, frequent cadence form. With the recent launch of SDO and continuing operation of SOLIS, such a database will be available to perform an analysis that could parallel the present paper. HMI produces 4096² pixel full-disk vector magnetograms as rapidly as every 90 s. In addition to allowing for a more detailed study of the photospheric magnetic field in space and time, this will open access to the transverse components where much of the enigmatic relationship to solar flares is likely to be found.

This work was supported by NASA MDI Grant NNX09AI90G. SOHO is a project of international cooperation between ESA and NASA. The authors thank the reviewer for insightful and helpful comments.