Following the same procedure as Garland et al. (2023) and Yurchyshyn et al. (2022), HMI magnetic field measurements were re-binned to 1 Mm pixel scale and transformed to a local Cartesian coordinate system using the same cylindrical equal area (CEA) projection used to produce the hmi.sharp_cea_720s series data (Bobra et al., 2014). The Cartesian boxes (or data cubes) in which the 3D magnetic fields would be extrapolated from the HMI magnetograms, which were the base maps of the boxes, were centered over the HARP center coordinates. In other words, the center of the base map was set to the HARP center coordinates. The box structures were quite large, spanning 400 by 200 by 200 pixels (re-binned) along the x-, y-, and z-axes (with z in the vertical) in order to accommodate as much of the field-line connectivity within the ARs and their surroundings as possible, which is a requirement for successful NLFFF application to real solar data (DeRosa et al., 2009). The specific size [400,200,200] accommodated the largest AR in the dataset. As described in the next section, specific portions of the field over the ARs within the data cubes were identified and used to calculate different parameters. Therefore, different data cube sizes for each flaring AR in the dataset did not need to be specified on a case-by-case basis, and a single size could be used to more easily automate parts of the modeling process so long as the field-line connectivity within the ARs was captured.

The HMI data retrieval, coordinate transformation, and NLFFF extrapolation were accomplished using the GX Simulator package (Nita et al., 2015; 2018; Nita et al., 2023). Following user specification of AR coordinates, observation time, and size/spatial resolution for the 3D data cube, the GX Simulator downloads the HMI magnetogram and remaps it to the bottom boundary of an empty data cube, and fills the empty cube with a potential field solution obtained from the normal component of the magnetic field at the lower boundary using a fast Fourier transform (FFT) method. In order to avoid unrealistic periodic, flux-balanced boundary conditions implied by the FFT solution at the lateral boundaries, more appropriate “open” boundaries are simulated by expanding the computational domain by 50% in each direction (Nita et al., 2023). The normal component of the field at the lower boundary is padded with a constant, non-zero value that is computed such that the unbalanced magnetic flux at the original lower boundary is compensated perfectly by the total signed flux from the added areas. The final potential field solution for the original domain is cut out from the expanded domain and used as the initial condition for the NLFFF extrapolation, which utilizes the weighted optimization code from Fleishman et al. (2017). During the extrapolation, the magnetic field at the top and side boundaries is frozen to be the potential field. Further details regarding the GX Simulator Automatic Model Production Pipeline (AMPP) can be found in the study by Nita et al. (2023).

For each of the 30 solar flare events in the dataset, the GX Simulator was used to generate [400,200,200] box structures containing the extrapolated magnetic field from the HMI magnetograms for a time period of 1 hour before through 1 hour after the flare start time. Given the 12-min cadence of the magnetograms, this meant that 11 3D solar magnetic field extrapolations were generated for each flare event. As mentioned previously, the box structures were centered over the HARP center coordinates, specifically, the coordinates from the HARP data time closest to the flare start time. The ±1 h window around each solar flare event was selected in order to specifically investigate the relative changes of the solar magnetic field shortly before and after eruption, when significant reconfiguration and energy release occur. Additionally, for the purposes of correlating the state of the field prior to the eruption and the different flare properties, this window of time limits the inclusion of data with an ongoing flare event happening ahead of the main flare being investigated. The ±1 h window also significantly limits the forecasting applicability compared to other studies that implement much longer preceding windows.

Similar to Garland et al. (2023), two regions were identified and used as masks to generate separate data volumes (i.e., specific portions of the total magnetic field within the box structures) from which the different magnetic field parameters would be calculated. The two regions were as follows:

1. Active region, AR

The definition of the AR in regards to the data volumes was the exact same as that used by Garland et al. (2023). The AR was defined as the areas where the magnitude of the observed photospheric vertical magnetic field (|B _z|) ≥ 300 Gauss (G) was within a specified bounding box of the magnetogram (see Figure 3). The bounding box was constructed by the following steps:

Step 1) Isolating pixels from the magnetogram with |B _z| ≥ 100 G.

Step 2) Running a smoothing mask over the filtered image.

Step 3) Identifying the bounded area that contains all pixels with

a pixel value ≥2 standard deviations above the mean.

For cases where multiple NOAA ARs are present within the FOV (such as in Figure 3), the bounding box is manually limited to the flaring AR. The threshold for |B _z| is slightly higher than the typical threshold range of 100–200 G seen in other studies to avoid noisy, weak field data (Bobra et al., 2014; Kazachenko et al., 2017; 2022). It should be mentioned that the threshold imposes the limitation that areas of highly twisted photospheric topology but weak flux, which provide information on the field line configuration emerging from the surface into the solar atmosphere (MacTaggart et al., 2021), are ignored. The implication of this is the assumption that all of the twisted topology critical to the flaring process are fully emerged at the time of flaring.

Using the above definition, the AR data volumes (or AR fields) were determined following the same method as Garland et al. (2023). Specifically, for each flare event, the AR at the HARP time closest to the flare start identified the photospheric footprints (or seed locations) from which the field lines within the [400,200,200] box were traced pointwise with a fourth order Runge–Kutta–Fehlberg algorithm, linearly interpolating between computational grid points, as part of the GX Simulator package (Nita et al., 2015; 2018; Nita et al., 2023). The coordinates of these field lines within the box specified the volume elements that would be used in the calculations of the magnetic field parameters; the collection of volume elements made up the data volume. The AR data volumes/fields were established (using the specified seed locations) for each of the 11 extrapolated magnetic fields throughout the 2-h window around each solar flare event.

The NPR replaced the overlap region (OR) from the study by Garland et al. (2023) and is defined as the region above the photosphere where the magnitude of the non-potential component of the magnetic field (|B _NP |) ≥ three standard deviations above the mean |B _NP | of the AR field, and either the unsigned twist number (|T _w|) ≥ one turn or the shear angle Ψ) ≥ 80°. The non-potential component of the NLFFF-extrapolated magnetic field was calculated as follows: B N P = B − B P o t , (1) where B is the NLFFF solution and B _Pot is the potential field solution. T _w is the twist number of the individual field lines calculated using the code developed by Liu et al. (2016), and the threshold of one turn is considered appropriate to define a moderate twist (Patsourakos et al., 2020). The shear angle was calculated as follows: Ψ = c o s − 1 B ⋅ B P o t | B ‖ B P o t | . (2) The 80° threshold for Ψ is the higher of two commonly used thresholds for shear (e.g., Lu et al., 1993; Leka and Barnes, 2007).

The purpose of the NPR was the same as that of the OR (Garland et al., 2023), which was to isolate the portions of the magnetic field expected to be involved with the flaring process and predict the location of the eruption. The hypothesis is that the portion of the solar magnetic field that has the greatest degree of non-potentiality above the photosphere would be associated with the location where energy is built up/released, that is, the location of the eruption of a solar flare. The thresholds for T_w and Ψ further stipulate that the non-potentiality of the field is due in large part to either or both of those properties being high, which is generally expected for flaring. The consideration of Ψ for the NPR, which was not considered for the OR, was motivated by the fact that both sheared magnetic arcades (SMAs) and magnetic flux ropes (MFRs), or a hybrid of the two, are the pre-eruptive magnetic field configurations considered viable (Patsourakos et al., 2020).

The NPR is established from similar parameters that are involved in the κ-scheme proposed by Kusano et al. (2020). The κ-scheme assumes a DA magnetic loop structure and provides a more specific estimate of flaring location (r _c < 1 Mm) along a PIL, whereas the NPR is a more generalized estimate (generally ∼ 10² Mm² in size) not constrained to a PIL and assumes the presence of a MFR, SMA, or hybrid configuration within the magnetic field structure.

Figure 4 shows a comparison of the OR and NPR for the X5.4 solar flare event that was analyzed as a case study by Garland et al. (2023). Both the OR and NPR aligned well with the area of initial flare brightening for this event. However, the NPR is less noisy and more concentrated over the point of eruption, whereas the OR contains many scattered, less coherent portions on the peripherals of the large sunspots. Figure 5 provides an example of an event where the NPR aligned far better with the area of flare brightening than the OR did, which was the common result for several test events from the full dataset. For the majority of the 30 flare events in the dataset, the NPR aligned well with the area of initial flare brightening (Figure 6). The alignment was worse for events where the AR did not have a well-defined PIL compared to those that did, as well as for events where the flare occurred on the edge of the AR away from an inversion line versus directly over one. In these cases, the magnetic field configurations were generally more conducive to emerging-flux-type flares (Thalmann, 2010) as opposed to that of “two-ribbon” flares. Details regarding how the areas of initial flare brightening were identified for the flare events using the SDO atmospheric imaging assembly (Lemen et al., 2012) 171 Å, 211 Å, and 1600 Å imagery can be found in the study by Garland et al. (2023).

Unlike the AR, which was defined at the surface and thus directly identified the seed locations that would be used to trace field lines from at each timestep, the NPR was defined above the photosphere within the corona. Therefore, the first step to determine the seed locations was to identify which field lines from the AR data volume/field came into contact with the NPR region at the start time of the flare. The footprints of those field lines determined the seed locations for the NPR data volume/field at each timestep throughout the 2-h window around the event (see Figure 7).

The magnetic field parameters calculated from the data volumes are shown in Table 2 with their description, units, and formulas. The summations for the parameters occur over the individual volume elements within the different data volumes, with the exception of ϕ, which is the surface flux. The parameters investigated were selected due to the potential predictive capability of these parameters found in other studies (e.g., Bobra and Couvidat, 2015; Kusano et al., 2020; Georgoulis et al., 2021). It is worth noting that all the investigated parameters are summed values, which was prompted by the findings of multiple studies showing that solar flares tend to be guided more by the total quantities than mean quantities (e.g., Leka and Barnes, 2007; Kazachenko et al., 2022). It is also worth highlighting that the parameters were calculated from data contained within the observed photospheric, modeled chromospheric, and modeled coronal fields. Whereas the finding of Korsós et al. (2020) shows that the optimal height for the WG_M method is 1,000–1,800 km, suggesting that the chromospheric magnetic field may be the key driver flare physics, signatures of flare occurrence have still been found within observed photospheric data and modeled coronal data. Therefore, data from all layers of the solar atmosphere were incorporated in the hopes of capturing the dynamics involved with the flaring process at each layer.

In addition to the parameters calculated directly from the different data volumes/fields themselves, the power spectra of the surface maps created for each of the parameters listed in Table 2 (with the exception of V) was calculated. This was motivated by the finding of Abramenko (2005) that the power spectrum of magnetograms is related to an AR’s future flare productivity/eruptive behavior. For each data volume/field (AR and NPR) at each timestep within the 2-h window around the flare events, surface maps were created for the magnetic field parameters by summing the volume elements along the vertical z axis. E _Free surface maps for both the AR and NPR fields for the X5.4 flare that occurred on 7 March 2012 are provided as an example in Figure 8. Following a similar approach as Abramenko (2005), the 2D power spectrum was calculated from each surface map. The total power was then calculated as the area under the power spectra curves within the wavenumber range of 0.3 rad/Mm ≤ k ≤ 3 rad/Mm, corresponding to the spatial scale 2Mm ≤ r ≤ 20 Mm, which is the ideal inertial subrange for magnetographic measurements (Abramenko, 2005). The NPR field E _Free power spectra for the X5.4 flare event are provided as an example in Figure 9.

An SPE plot is a sensitive tool used to discover the average (possibly weak) response of a complex system to a particular event (Mason and Hoeksema, 2010). The SPE plots are combinations of the individual time series for the 30 flare events in the dataset, aligned by the flare start times. Unlike the time series plotted for the case study analysis (Garland et al., 2023), which were generated from the parameters normalized to the value at the flare start time, the individual time series used in the SPE plots were generated from the field parameters normalized to the value at the start of the 2-h window around the flare event. Such normalization still allowed for effective comparison of the relative changes, but it also revealed the consistency in any trends from a forecasting perspective. The data shown in the plots are the mean and standard deviation of the individual time series.

SPE plots were produced for each data volume/field for the solar flare dataset as a whole, as well as for flare events of similar magnetic and sunspot classifications. However, as Table 1 indicates, the categorization by each specific magnetic and sunspot classifications would result in too small data subsets. In order to maintain a degree of statistical significance, data subsets were established for only the most common Mount Wilson magnetic classification and modified Zürich/McIntosh sunspot classification. Thus, SPE plots were produced for the following:

1. Full solar flare dataset—30 flare events (18 X-class and 12 M-class)

The SPE plots for the full flare dataset are shown in Figures 10 and 11. Overall, the parameters that showed the most consistent trends among the flare events, particularly for the hour leading up to the eruption, were ϕ, E _Free , τ _Tot, and ρ _Tot. Consistency was based on the relative standard deviation, which was calculated as the standard deviation divided by the range of the averaged data. The parameters ϕ, E _Free , τ _Tot, and ρ _Tot had the four lowest relative standard deviations, respectively, and therefore will be the focus of discussion for this section. The trends seen for these parameters were the same across all the data subsets. However, it is worth noting that a greater relative change and smaller standard deviations were found for the subsets of flare events with any magnetic classification other than βγδ and any penumbra/distribution classification other than ‘kc’. Although this could just be a result of the smaller sample sizes, it could also be due to the nature of the larger and more complex ARs, which showed greater overall variability in the evolution of the parameters, particularly after eruption. In several of the flare events with the largest and most complex ARs, the post-flare values of the parameters quickly increased back to the pre-flare values, indicating that the magnetic fields did not fully relax to less stressed states. These ARs were typically more active, in line with the well-known rule of thumb of flare persistence (Mason and Hoeksema, 2010).

The most consistent trend across all the flare events was the steady increase in ϕ for the AR field over the hour leading up to the flare start, peaking at the time of the flare, and then decreasing throughout the remainder of the event. This could be explained by the processes of flux emergence and cancellation, which have been associated with flaring activity (e.g., Nitta et al., 1996; Li et al., 2000; Wang et al., 2002; Kutsenko et al., 2021). It could also be the result of the stretching (in the vertical) and relaxing of the magnetic field for the “two-ribbon” flare model. As the strapping or envelope field lying over the expanding flux rope is stretched in the vertical, the z-component of the field and therefore ϕ increases. Once the instability of the magnetic field reaches the point to trigger the eruption, the magnetic field relaxes to a less stressed state in the form of the simple post-flare arcade that is no longer stretched as extensively. Leka and Barnes (2003) studied the temporal evolution of the ϕ during a similar time interval for two ARs that produced X- and M-class flares, finding “some flux growth” prior to the flares, but no post-event decrease. The difference in the observed trend for ϕ found in this study and theirs could be due to the different magnetogram data used or the higher AR defining threshold (|B _z|≥ 300 G) utilized in the present study. Similar to that the values in the study by Leka and Barnes (2003), the relative changes in ϕ are small, only about 1%–5%.

An opposite trend is seen for ϕ for the NPR field, where on average, there was approximately 10% decrease over the two timesteps leading up to the start of a flare, followed by a roughly equal increase over one to two timesteps after the flare start (Figure 11). Relating to the changing structure of the magnetic field, this can be explained by the shearing and reconnecting of the non-potential field lines associated with the SMA/MFR structures as the field destabilizes, leading up to the flare, and then relaxes into the post-flare arcade after the eruption. Essentially, field lines in close proximity to the expanding flux rope reconnect below the flux rope, causing them to wrap around the twisted axis. As these field lines become a part of the flux rope, they become more elongated in the horizontal, leading to a decrease in the z-component of that portion of the magnetic field. Following the eruption and decay of the flux rope, the post-flare arcade sets up with more vertically oriented, simple loop structures compared to the flux rope that was present prior to eruption. Higher cadence data would be useful to further examine the rate of this “dip” in ϕ for the NPR field around the start of the flare.

On average, for the NPR field, there is an initial increase in E _Free until 24 min before flare start, followed by a slight decrease up to the flare start and a more significant decrease from the flare start to 1 h after. The mean dissipated energy is 41% ± 19% of the maximum free energy, which is lower yet overlapping with the 60% ± 26% found by Aschwanden (2019). The initial decrease in E _Free following the peak just before the flare start suggests that the free magnetic pressure reaches a point where the field can no longer fully contain it and it starts to “leak” out. This “leaking” of the pressure and slight loss of free energy is likely brought about by the ongoing reconnection as the field destabilizes but before the critical point of destabilization is reached and the bulk of the energy is released. In other words, once enough free energy has been built up and the free or non-potential magnetic pressure is high enough, each reconnection event that occurs before the eruptive instability is reached releases a small amount of the free energy. It is plausible that these small releases of the free energy just before the eruption are associated with the short-lived small-scale brightening seen in UV and EUV imagery (e.g., Leka et al., 2023).

It should be noted that while the trend just discussed was the average trend seen in many of the 30 flare events investigated, some events did show more a more variable evolution of E _Free , with intermittent instances of energy buildup and dissipation superimposed on the overall change of higher pre-flare energy to lower post-flare energy. Other studies that have analyzed the temporal evolution of E _Free estimated from NLFFF extrapolation have reported only step function decreases in E _Free following the occurrence of the strong flares examined (Thalmann and Wiegelmann, 2008; Murray et al., 2013; Sun et al., 2012), in line with the store-and-release concept for solar eruptions. However, the findings of Aschwanden (2019) and the present study suggest that the evolution of energy during solar flares is more complex than a single-step function, exhibiting a more fragmentary and impulsive nature. The different findings could be the result of different temporal resolution, extrapolation techniques, and/or the portion of the modeled 3D solar magnetic field from which E _Free was calculated.

The evolution of τ _Tot and ρ _Tot for the NPR field, on average, is similar to that of E _Free , which was expected. For the majority of flare events examined, |T _w| either increased or remained fairly constant in the hour before flare start and then decreased following the eruption as the field relaxed. This was more evident in the flare events where highly twisted (|T _w| ≥ 2) structures were present in the pre-flare field configuration than in events where the pre-flare twist was more moderate (|T _w| ≈ 1), which saw a more variable evolution of twist. Considering τ _Tot and ρ _Tot are essentially T w Tot weighted by |B _z| and |B _NP|, respectively, the evolution of τ _Tot and ρ _Tot is also a factor in the evolution of the vertical and non-potential fields, which is indicated by E B z and E _Free .

Similar to the study by Kazachenko et al. (2017), the relationship between different flare properties and the magnetic field parameters was quantitatively described by the Spearman ranking correlation coefficient (r _s ), which measures the monotonic relationship between two variables. The different flare properties are as follows:

1. GOES measured peak X-ray flux in W/m ².

The method of bootstrapping (Efron and Tibshirani, 1986) was used to estimate the uncertainty in the correlation coefficient due to sampling. Specifically, r _s was calculated 10⁴ times for 30 randomly selected data pairs from the dataset with replacement. The mean r _s calculated from the bootstrap method was the reported r _s for the particular parameter and flare property correlation, with the standard deviation of r _s reported as the correlation’s statistical uncertainty. In addition to r _s , the p-value (p) and its uncertainty were calculated as well. The commonly used thresholds for r _s to express the qualitative strength of the correlation are as follows: |r _s | ∈ [0.0, 0.4)—weak, |r _s | ∈ [0.4, 0.6)—moderate, |r _s | ∈ [0.6, 0.8)—strong, and |r _s | ∈ [0.8, 1.0]—very strong (Dodge, 2008). For strong and very strong correlation coefficients, the relationship between the parameter and flare property was fitted with a power-law function of the form Y = a X b , (3) where Y is the flare property and X is the magnetic field parameter. The scaling factor, a, and exponent, b, were found using linear regression in the log–log space.

The process of correlation coefficient calculation was carried out for each data volume/field (AR and NPR), each set of parameters (field parameters and power spectrum parameters), and each timestep leading up to flare start, as well as the sum across the timesteps leading to the start. Tables 3 and 4 provide a summary of r _s values calculated from the AR data volumes for the field parameters and power spectrum parameters, respectively. Tables 5 and 6 provide the same for the NPR data volumes. Listed in the tables are the maximum r _s values from the different timesteps for each parameter and flare property examined.

In this study, the 3D solar magnetic field was modeled using optimization-based NLFFF extrapolation during 30 strong and extreme solar flare events. Specifically, 3D magnetic fields were modeled using the 12-min cadence HMI photospheric vector magnetograms, spanning a time period of 1 hour before through 1 hour after the start of 18 X-class flares and 12 M-class (greater than M6.5) flares. Using the modeled data, several magnetic field parameters were calculated from the fields directly, as well as from the power spectrum of surface maps generated by summing the fields along the vertical axis, for the following two different regions:

1. Active region (AR): areas with photospheric |B _z|≥ 300 G.

The superposed epoch (SPE) plots of the magnetic field parameters were analyzed to investigate the evolution of the 3D solar field during the solar flare event and discern consistent trends across all solar flare events in the dataset, as well as across subsets of flare events categorized by their Mount Wilson magnetic classification and modified Zürich/McIntosh sunspot classification. The relationship between the different flare properties (peak X-ray flux, duration, rise time, decay time, impulsiveness, and integrated flux) and the magnetic field parameters was quantitatively described by the Spearman ranking correlation coefficient, r _s . In addition to the individual parameters, ratios of the parameters to the free energy, as well as the net current helicity of the AR fields weighted by those ratios, were included in the list of parameters from which correlation with the flare properties were calculated. The key results of the study are as follows:

1. For the majority of the 30 flare events in the dataset, the NPR aligned well with the area of initial flare brightening through the start of the eruptive/flash phase. The alignment was worse for events where the magnetic field configuration was more conducive to emerging-flux type flares, that is, the AR did not have a well-defined PIL or the flare occurred on the edge of the AR away from an inversion line versus directly over one. Overall, the NPR provides a good indication of where within the magnetic field over an AR an eruption will occur if triggered.

The findings of this study suggest that the specific magnetic field parameters calculated from the identified portions of the modeled 3D magnetic field have potential predictive capability and could benefit flare forecasting efforts. The location of a flare within an AR if triggered can be forecasted reasonably well through identifying the NPR. Although a perfectly consistent, distinct signature that specifies the exact timing of flare occurrence may not be seen in the time series of the parameters calculated from the modeled magnetic field, the destabilization of the field that is conducive to an eruption is observable. The relationships discussed in this study, calculated solely from data that were within an hour timeframe leading up to flare occurrence, could provide reasonable estimates of the flare properties of any strong/extreme flare that may occur within an hour’s time.

While the trends and relationships found from the statistical analysis show high potential for flare forecasting, a major limitation of the study is the number of events included in the dataset. Although as a rule of thumb a minimum of 30 observations is generally considered sufficient to conduct significant statistics, the confidence in the results is limited. Thus, the next logical step would be to expand the dataset, initially including more X- and M-class flares, then expanding into the lower-class flares as well to find possible differences between the weaker and stronger flares. Following the initial expansion of more flare events, the next step would be to include data from non-flaring ARs to test whether the trends that hold consistent across flare events are discriminating in regards to flare occurrence.

Another limitation of the study is the use of only a single model; therefore, it was not able to determine if the findings are an artifact of the specific model used. To make such a determination, the analysis could be re-performed following the exact same procedures, only utilizing different models to extrapolate the 3D solar magnetic field. If the trends and relationships found continue to hold for varying extrapolation techniques, it would increase the confidence that the findings are grounded in the actual physical processes occurring within the solar magnetic field during flare events. Of particular interest would be the use of models that utilize physics-informed neural networks, such as the method of Jarolim et al. (2024), which builds upon the work of Jarolim et al. (2023) by incorporating chromospheric magnetic field measurements into the extrapolation process to generate more realistic field simulations. The inclusion of the chromospheric data helps correct the lower portions of the solar magnetic field, where the force-free assumption is not valid due to the dense plasma conditions present. As discussed in the study by Korsós et al. (2024), such neural network-based methods provide more realistic field extrapolations with considerable improvement in computational efficiency compared to the standard NLFFF extrapolation, proving more suitable with respect to real-time flare forecasting application.

Finally, the 12-min cadence between observations of the standard hmi.B_720s dataset brings about another limitation: for time scales less than 12 min, for example, flare eruptions, changes in the magnetic field are lost. Sun et al. (2017) created a new high-cadence vector magnetogram dataset for HMI that is more suited for investigating the evolution of the magnetic field during such quick phenomena. The 135-s cadence series hmi.B_135s has the same format and is processed with nearly identical pipeline options as the 720-s data (Sun et al., 2017). Although the noise level in the field strength for the higher cadence dataset is about 50 G higher than that of the standard dataset, measurements in the strong-field regions (|B|≥ 300 G) agree well between the two. Performing the same analysis with the higher cadence could highlight signatures that were potentially missed with the standard 12-min cadence dataset.