© C. Guennou et al., Published by EDP Sciences 2017

1. Introduction

Solar flares are sudden brightenings of the solar atmosphere occurring over the whole spectral range, from radio to X-rays, resulting from the release of a huge amount of magnetic energy. Flares can be classified in terms of the strength of the X-ray flux measured close to Earth. They can be divided into two general subcategories: eruptive flares, accompanied by a coronal mass ejection (CME), large ejection of magnetized plasma released into the interplanetary space, and confined flares, for which only electromagnetic radiation is emitted. Both CMEs and solar flares can impact the near-Earth environment and in particular human technologies. The most energetic solar flares can heat and ionize the upper Earth atmosphere, engendering disruption in radio communication and Global Positioning System (GPS) inaccuracies. CMEs can be accompanied by the burst of energetic particles and can be associated with strong geomagnetic storms, inducing ground-level electric fields likely to produce power grid damages. Astronauts, pilots, and satellites can also be affected by radiations originating from flares or CMEs. The ability to accurately predict both solar flares and CMEs is therefore fundamental to protect both space and ground-based technologies from strong space-weather events.

The European H2020 research project FLARECAST (Flare Likelihood and Region Eruption Forecasting – http://flarecast.eu/) aims to develop a fully automated solar flare forecasting system. FLARECAST will automatically extract various magnetic field parameters of solar active regions (ARs) from vector solar magnetograms to produce accurate predictions using the state-of-the-art forecasting techniques based on data-mining and machine learning. Various systems of prediction invoking different categories of models have been developed in the past decades, as e.g., the “Theophrastus” tool (McIntosh 1990), the linear-prediction system of Gallagher et al. (2002) used by SolarMonitor, the discriminant analysis of Barnes et al. (2007), the Automated Solar Activity Prediction (ASAP) of Colak & Qahwaji (2009), based on machine learning, and more recently the statistical learning technique of Yuan et al. (2010). A recent comparison between the current forecasting tools using line-of-sight (LOS) magnetograms has been performed by Barnes et al. (2016), showing that none of them substantially outperformed all others. All these forecasting techniques, including the future FLARECAST forecasting tool, require scalar quantities derived from photospheric magnetic field to be able to make flare and/or eruption predictions. In this paper, we focus on the predictability capabilities of these scalar quantities, in order to determine the most reliable for flare and eruption forecasting use.

It is now widely accepted that the energy source of solar flares is stored in highly non-potential magnetic fields. The storage of energy, typically estimated in the 10²⁸–10³² erg range (Schrijver et al. 2012), arises from a long phase of magnetic stress and free magnetic energy build-up before sudden flare or/and eruption. A large variety of models have been recently developed to explain this storage and released mechanism (see e.g., the recent reviews of Aulanier 2014; Janvier et al. 2015; Lin et al. 2015; Schmieder et al. 2015), but it is still unclear why some ARs trigger a flare, while others remain quiet. In this context of absence of a clear physical scenario for flare triggering, efforts have been concentrated on empirical flare prediction methods, investigating the behavioral patterns of ARs.

Mapping the 3D coronal magnetic field is not systematically practicable, while the photospheric magnetic field is more easily measurable using spectro-polarimeters such as the Helioseismic and Magnetic Imager (HMI; Schou et al. 2012) on board the Solar Dynamics Observatory (SDO). Consequently, flare forecasting tools are mainly based on quantities derived directly from solar magnetograms. As flare occurrence is correlated with size, variability, and complexity of the ARs, i.e., the associated non-potentiality degree, many attempts have been made in order to find reliable photospheric eruptivity indicators. Global photospheric features, such as magnetic field gradient, currents, magnetic geometry, and magnetic free energy, have been proposed to establish a link between photospheric observations and coronal activity.

However, it has been shown that observing clear pre-signature of flare and eruption in observations is challenging. Many studies investigated the link between changes in some photospheric parameters and the coronal eruptive and flaring activity, by investigating prior temporal changes, performing superposed epoch analysis, or forecasting the likelihood of the flaring events (see e.g., Bao et al. 1999; Leka & Barnes, 2003a; Schrijver, 2007; Jing et al. 2010; Mason & Hoeksema, 2010; Falconer et al. 2011; Ahmed et al. 2013; Bobra & Ilonidis, 2016, and references therein), with moderate results. More recently, Al-Ghraibah et al. (2015) studied the magnetic parameters of about 2000 ARs to search for the best eruptive predictor. They used a wide range of parameters, including a wavelet analysis to resolve multiple-scale changes, and found that that magnetic field properties alone are not sufficient for powerful flare forecasting, although the magnetic-gradient related features appear to be one of the most reliable indicators for flare predictions.

Despite all these previous studies, none of the current photospheric indicators, or combinations thereof, have yet demonstrated an unambiguous eruptive or flaring criterion. However, controlled cases (e.g., originating from numerical datasets) have barely been used to test the predictability of these parameters. Kusano et al. (2012) presented a parametric analysis of eruption onset, by varying two parameters associated with the magnetic structure. They showed that these two parameters are able to discriminate between eruptive and non-eruptive ARs, although no comprehensive scalar quantity directly measurable is provided. In this work, we use MHD numerical simulations of the formation of stable and unstable magnetic flux ropes (Leake et al. 2013, 2014) in order to systematically investigate the pre-eruptive signature included in different magnetic parameters. This series of numerical experiments is based on the emergence of a convection zone magnetic flux tube into a solar atmosphere. The interaction of the emerging magnetic flux with the pre-existing magnetic field plays a key role for triggering impulsive events. This class of simulations is thus able to explain and reproduce multiple active solar phenomena (see e.g., the review of Cheung & Isobe 2014). Time series of magnetograms from parametric simulations of stable and unstable flux emergence, i.e., corresponding to respectively quiet AR versus eruptive flare configurations, are used to compute a large range of parameters. This list includes parameters previously used for operational forecasting, as well as parameters used for the first time for this purpose, such as helicity (Pariat et al. 2017), the current-weighted magnetic inversion line, and the length of the strong-current magnetic inversion line portions.

The remainder of this paper is as follows. The set of eruptive and non-eruptive simulations is summarized in Section 2. Section 3 describes the magnetograms and the magnetic parameters extracted from our simulation sets. The results are presented in Section 4, and a discussion on the impact of data masking on the analysis is proposed in Section 5. A basic analysis of the noise impact on our results is presented in Section 6. Section 7 provides a parametric study of the inversion line related parameters and the influence of the different thresholds on eruption predictability. Finally, Section 8 summarizes our work and presents our main conclusions.

2. Numerical datasets: parametric simulations of eruptive and non-eruptive active regions

2.1. MHD simulations

In order to test the reliability of photospheric eruptive predictors, we used the three-dimensional visco-resistive magnetohydrodynamic (MHD) simulations of Leake et al. (2013, 2014). These simulations describe the partial emergence process of a twisted magnetic flux tube into a stratified solar atmosphere, where a coronal arcade field is present. Both stable and non-stable flux ropes are formed as a result of the flux emergence, depending on the choice of certain parameters. These simulations are thus analogous to respectively quiet ARs and eruptive flare productive ARs, where the newly formed flux ropes are ejected higher in the solar corona.

The evolution of the system is described by the visco-resistive MHD equations (see Eqs. (1)–(4) from Leake et al. 2013), and the plasma is assumed to be fully ionized. The MHD equations are solved using the Lagrangian-remap code Lare3D (Arber et al. 2001), using an irregular Cartesian grid. The initial conditions consist of a hydrostatic background atmosphere, stratified such as the solar convection zone, the photosphere/chromosphere, the transition region, and the corona. An arcade field covering the entire simulation domain is imposed on this background atmosphere (cyan solid lines in Fig. 1), and a right-hand twisted flux tube is inserted in the solar convection zone (red solid lines in Fig. 1), aligned along the y axis. The arcade field is transitionally invariant along the y axis, generated by a source much deeper in the solar convection zone than the initial horizontal flux tube. This background coronal field is designed to reproduce the magnetic field of an old decaying active region. In these simulations, the initially buoyant twisted flux tube, which is line-tied at the side boundaries, partially emerges from the convection zone into the corona due to a perturbation in the flux tube’s pressure and density (see Eqs. 18 and 19 from Leake et al. 2013) at the center of the flux tube. The later evolution of the emerging field involves surface shearing and rotation, and the formation of a new coronal flux rope above the new active region. The ultimate state of this new flux rope depends on the choice of the overlying field parameters.

Fig. 1. Emergence of the convection magnetic flux tube into the solar atmosphere at t = 110 t 0 for both the stable SD (A) and the unstable SD E (B) simulations. The cyan lines trace the coronal arcade magnetic field, originating from the lower boundaries, whereas the red lines belong to the magnetic flux tube and originate from the y = ±maxy side boundaries. The gray scale indicates the magnetic field strength at the surface.

The emergence process of the flux tube is examined through a range of initial coronal arcade strengths and orientations, in order to investigate the effects of the external magnetic field on the formation mechanisms of both stable and unstable flux ropes. The flux rope, self-consistently formed by the emergence of the flux tube, is either ejected, i.e., corresponding to an eruptive flare, or confined in the corona, i.e., equivalent to a quiet stable AR, depending on the orientation of the coronal arcade. For the non-eruptive simulations, the direction of the arcade field is the same as that of the top of the flux tube, maximizing the confinement of the flux tube by the arcade field, as shown in Figure 1A. Conversely, the eruptive simulations are driven by the magnetic reconnection between the arcade and the top of the flux tube’s axis, due to the opposite orientation of the magnetic field in the two structures. In this case, the horizontal magnetic field B_x changes sign at the interface between the flux tube and the arcade, making this separatrix a favorable location for magnetic reconnection. For each of the two arcade orientations, the arcade magnetic field strength is varied, leading to three different simulations: strong (SD), medium (MD), and weak (WD), corresponding to the initial surface field strength of 26, 19.5, and 13 G. An additional simulation where no coronal field is present (ND) is also performed, resulting in the formation of a stable flux rope.

For the eruptive simulations, the coronal magnetic field strength varies in the same way, namely {WD E, MD E, SD E}, where E corresponds to “Eruptive” simulations. The continued reconnection process between the emerging field and the arcade starts as the flux rope reaches the corona, changing the connectivity of the system and forming magnetic field lobes on both sides of the emerging flux rope. Figure 1B shows the magnetic field lines configuration of the flux tube emergence process for the SD E simulation, at t = 110 t₀ (with t₀ = 55.7 s, see Sect. 2.2 for details), slightly before the ejection of the flux rope. The formation of a shearing quadrupole configuration above the photosphere is clearly visible, with a central arcade of the emerging structure (red lines) now allowed to expand higher in the corona, whereas horizontal expansion is limited by the lobes on either sides (in red). Once the central arcade reaches a sufficient height, internal reconnection takes place at about t ~ 120 t₀ (depending on the overlying arcade strength, see Leake et al. 2014, for details) beneath the flux rope axis, and the newly formed flux rope accelerates its vertical expansion, and the flux rope is immediately ejected. The objective of this work is to investigate if some pre-eruptive flare variations in one or more physical parameters are detectable using only photospheric magnetograms, at a reasonable stage before t ~ 120 t₀, providing thus a reliable eruptive indicator.

This set of seven simulations treats the cases of both eruptive and non-eruptive flux rope formation. The presence and the orientation of the coronal arcade are critical parameters for the eruption onset. The ratio of the arcade flux to emerging flux is also a fundamental parameter, controlling the flux rope vertical acceleration, its size, and the amount of reconnection allowed to occur.

2.2. Scaling

The visco-resistive MHD equations are non-dimensionalized, using a normalization constant for each variable. In this work, we choose a slightly different normalization than the original paper, in order to increase the representative size and the magnetic flux of the active regions, to obtain more representative ARs. For the length, we rescale the simulations using the normalizing value L₀ = 8.5 × 10⁵ m, corresponding to five times the original value adopted by Leake et al. (2013, 2014; see Sect. 2.2). The two other normalizing values, the magnetic field (B₀ = 1300 G) and gravitational acceleration (g₀ = g_sun = 274 ms⁻²), remain unchanged. The other normalizing variables that are affected by this rescaling are, namely the density (ρ₀ = 5.77 × 10⁻⁵ kg m⁻³), the velocity (v₀ = 15.26 × 10³ ms⁻¹), the time (t₀ = 55.7 s), the temperature (T₀ = 28.2 × 10³ K), the current density (j₀ = 0.122 Am⁻²), the viscosity (ν₀ = 7.49 × 10⁵ kg m⁻¹s⁻¹), and the resistivity (η₀ = 16.3 × 10³ Ωm). The Reynold’s number and the magnetic Reynolds number are kept in the same order as in the initial simulations, both evaluated to 100.

Using a larger spatial scaling than that of the initial analysis allows us to obtain larger active regions and larger magnetic flux, therefore more representative of the observed eruptive active region. Still, we are limited to the study of a relatively small active region with a characteristic flux of 1–2 × 10²¹ Mx during the early emergence process, and a characteristic size of about 30 Mm, roughly giving an area of about 296 μhs (i.e., microhemisphere – 1 μhs = 3.04 × 10⁶ km²). According to Sammis et al. (2000), our simulations belong to the smallest flaring active region. For comparison, the authors have shown that during the 1989–1997 time frame, all the X4 flares only originated from active regions with an area greater than 1000 μhs. It should be noted that not all active regions in this size range are flaring: size seems to be a necessary condition to X4 flares, but not sufficient. Given our framework, our study is thus limited in terms of size and complexity of active regions that we cannot explore: the reliability of eruptive indicators is tested only for a given class of active region size and complexity. Nonetheless, it is worth noting that a recent study by Toriumi et al. (2017) concludes that a complex δ-sunspot is not a necessary condition for flaring ARs, even for the X-class flares. Thus, this kind of controlled study is important in the sense that the size and the complexity of an active region are not discriminant parameters for eruptivity: small active regions can still produce flares, while two active regions in the same size and complexity range do not necessarily have the same flaring likelihood.

However, if the rescaling of the simulation set enables us to have greater size, it also influences other quantities used in the simulations. Using our scaling, the transition region is 8.5 Mm thick (1.7 Mm in the original simulations), which is significantly higher than its characteristic size of about 0.1 Mm. However, given the spatial resolution imposed by the total domain size and the computational cost, such a thickness is required to resolve the large temperature and density gradient occurring in this part of the solar atmosphere.

3. Analysis: extracting the physical properties of ARs

In the following, the methodology that we follow to investigate the reliability of the eruptive flare indicators is detailed. From the series of 2D-plane vector magnetograms of the three eruptive and four stable simulations, we compute a series of parameters in order to detect some particular behavior associated only with the eruptive simulations. A predictor can be considered as reliable if it shows significant change(s) in a reasonable timescale prior to the eruption. Since flares and eruptions are driven by a store-and-release mechanism of magnetic energy, eruptive indicators must provide indications of such process occurring. Such signature includes threshold beyond which the system becomes unstable, most likely generating an eruption or a flare. Therefore a reliable predictor should present a different behavior respective to the eruptive or the non-eruptive nature of the simulations, but a similar behavior for simulations of the same nature. This signature should also be significant enough to be measured, and detectable in a sufficient time prior to the eruption to be able to then perform operational forecasting. Here, we focus only on determining which quantities could potentially be proficient eruptive flare predictors, i.e., parameters associated with higher values for the three eruptive simulations {WD E, MD E, SD E} than that of the stable {ND, WD, MD, SD} prior to the eruption. The associated thresholds needed for performing operational predictions require the use of (1) real photospheric magnetogram observations and (2) substantially more AR samples to be statistically significant.

3.1. Magnetograms

From the time series of the 3D MHD simulations, we first extract photospheric-like 2D-plane vector magnetograms, by interpolating the cube domain at z = 0 (as defined by Leake et al. 2013), onto a regular grid using a resolution of 0.86 L₀ (730 km pixel⁻¹) in both directions, which is equivalent to the instrument resolution of about 1′′. This resolution is the same as that of the SDO/HMI magnetograms. In our simulations, the surface (z = 0) is defined as the beginning of the temperature minimum region, i.e., corresponding to the region of the lowest pressure scale height. Since the eruptive flare occurs around t ~ 120 t₀ for the three eruptive simulations {WD E, MD E, SD E}, we restrict the time windows of our analysis from t = 0 t₀ to t = 150 t₀, using a time sampling of Δt = 5 t₀, where t₀ = 55.7 s. Therefore, we finally obtain 31 magnetogram series for each vector magnetic field component B_x, B_y, B_z, and for each of the seven parametric numerical simulations {ND, WD, WD E, MD, MD E, SD, SD E}.

We first apply a mask to the entire magnetogram series, excluding the pixels for which , where B_mask = 30 G. This B_mask threshold has been chosen in order to be greater than the background magnetic field due to the coronal arcade field, respectively, about 26, 19.5, and 13 G for the SD, MD, and WD simulations (see Sect. 2.1). In observational analysis, a threshold is usually applied to exclude noisy data and only retain strong field areas. Our data mask threshold is relatively low compared to the typical uncertainties measured in observed magnetograms and so far the usual threshold applied to data. This is due to the relatively small typical values of magnetic field in our simulations, with B_z values in the [−800 G, 800 G] range, due to the small size and flux of our ARs (see Sect. 2.2). Many different masking methods have been used in previous analyses, using different threshold values applied on either B or B_z magnetograms, and will be briefly introduced in Section 5. The impact of applying a data mask on the detection of photospheric eruptive signatures will also be discussed in more detail in Section 5. The effects of noise, using a random perturbation of the magnetograms, will also be investigated in Section 6.

Figure 2 shows an example of such magnetograms, for both eruptive and non-eruptive simulations. Figure 2A displays B_z magnetograms for the MD simulation, either masking (right) or not (left) of the data, while Figure 2B exhibits the same photospheric maps for the MD E simulations. Masking data has a great impact on the area of the region of interest, used then for the computation of the different eruptive indicators. Eruptive and non-eruptive magnetograms are very similar, except near the external edge of the polarities. This is due to the quadrupolar nature of the eruptive simulations: as mentioned before, the coronal arcade has an opposite orientation with respect to the flux tube, i.e., the arcade B_x component has the opposite sign, to favor reconnection above the flux tube axis. From such time series of masked magnetograms, we compute a set of scalar parameters, most of them typically used in standard solar flare operational forecasting methods.

Fig. 2. Magnetograms of the Bz component at t = 100 t0 for both the MD (A) and MD E (B) simulations. Both non-masking (left) and masking (right) magnetograms are displayed, in order to show the reduction of the region of interest for further analysis. Eruptive and non-eruptive simulations present very similar magnetograms, except near the external edge of the polarities, due to the quadrupolar geometry above the surface.

3.2. The parameters

All the quantities considered in this work are scalar and can be derived using exclusively vector magnetic field data. They described the overall physical condition of an active region and the ongoing evolution of its magnetic field. Most of them have been chosen based on their previous identification as potential flare predictors, but new quantities such as helicity or WL_sc, the current integral along the magnetic polarity inversion line, are also tested. Each derived quantity is parametrized such as one single number is able to characterize the state of the whole active region at a given time. For quantities that are spatially distributed, such as B_z (x, y), we compute the four first moments as described in Leka & Barnes (2003a): the mean, the standard deviation, the skewness, and the excess kurtosis (i.e., the kurtosis −3, for comparison to normal distribution). The skewness describes the asymmetry of the spatial distribution, while the kurtosis accounts for small changes in extremal values. As pointed out by Leka & Barnes (2003a), computing the four first moments of the spatial distribution allows us to quantify the subtle evolution of the different parameters.

3.2.1. Magnetic field

This first class of parameters is directly related to magnetic field evolution. All these parameters, listed in Table 1, come from the work of Leka & Barnes (2003a, hereafter LB03), and the readers are referred to Section 3 of their paper for a full detailed description. The spatial/scalar column specifies if the computed parameters are spatially distributed or not. If so, the four first moments are computed, easily identifiable by the notation. The temporal evolution of the total, vertical, and horizontal magnetic fields, respectively, noted B, B_z, B_h, and their associated moments provide insight about how the emerging field changes the distribution of B, and the overall polarity imbalance of the active region. The associated horizontal spatial gradients , and the corresponding moments indicate how much the magnetic field is sheared and distorted, an important indicator of the AR non-potentiality, especially along the neutral line. We also compute four quantities related to the AR flux: the total unsigned flux ϕ_tot, the net flux ϕ_net, and the positive and negative flux Φ_±, providing insights into the overall magnetic flux and the flux imbalance during the partial emergence process. Finally, we compute the photospheric free magnetic energy, which refers to the energy directly available for eruptive activity, i.e., measured with respect to the potential field energy (see the corresponding equation in Table 1). It is worth noting that in our case, the free magnetic energy is computed from surface integrals, not from the whole AR volume, which would require 3D magnetic field extrapolations.

3.2.2. Magnetic field geometry

This second class of parameters, listed in Table 2, quantifies the morphology of the AR magnetic field and its deviation from the potential configuration. The inclination angle γ measures how much a magnetic field is inclined relative to the vertical axis. It reflects the magnetic field orientation, providing hints about the flux emergence evolution. In the present work, we adopt the following definition(1)

As defined, the inclination angle tends toward small values when the magnetic field is approaching the vertical, while great values correspond to an almost horizontal magnetic field. This quantity indicates the morphology evolution of the magnetic field during the emerging process.

In a force-free field framework, the Lorentz force is null and the magnetic field follows(2)where α is the local twist parameter, referring to the torsion of each individual field line. Using magnetogram data, α can be approximated by (see e.g., Pevtsov et al. 1994)(3)where μ₀ is the permeability in free space. Using a set of 133 of flaring ARs, Nindos & Andrews (2004) found that the pre-flare value of the twist parameter α was, in general, higher for M-class eruptive flares, suggesting the twist parameter as a reliable indicator for CMEs.

The magnetic shear angle has been extensively used in previous studies to measure the non-potentiality of active regions. Initially introduced by Lu et al. (1993) and Hagyard et al. (1984), the three-dimensional shear angle Ψ measures the angle between the observed magnetic field and its potential component, while its horizontal projection, namely the “horizontal” or “planar” shear angle , quantifies the azimuthal difference between the observed and potential magnetic field. The potential magnetic field is computed following the method of (Valori et al. 2012). We also tested the additional shear angle-related parameters A[Ψ > 80°] and A[ > 80°] proposed by Leka & Barnes (2003a), corresponding to the AR area where respectively the 3D and the horizontal shear angle exceeds 80°.

Magnetic helicity measures to what extent the magnetic field lines are wrapped around each other, and how much the individual magnetic field lines are twisted and writhed, relatively to their lowest energy state. This parameter provides a quantitative estimation of the geometric properties of the magnetic field lines. Because helicity is a conserved MHD quantity, even in resistive MHD where the dissipation is very small (Pariat et al. 2015), im portant efforts have been carried out concerning its estimation and its relation with solar flares (see e.g., Nindos & Andrews 2004; Démoulin 2007; Démoulin & Pariat 2009; Park et al. 2010). In the present study, helicity has been computed using the G_θ proxy for the helicity flux density from Pariat et al. (2005), using the series of photospheric magnetograms. In this framework, the time variation of the relative magnetic helicity can be written as(4)where is the photospheric surface, B_z is the magnetic field normal to this surface, r = x – x′ is the vector between the two photospheric positions x and x′, with u and u′ the associated flux transport velocities of the photospheric field line footpoints. While the flux transport velocity could be directly extracted from the simulations, the velocity u has been estimated using the differential affine velocity estimator for vector magnetograms (DAVE4VM) from Schuck (2008) in order to analyze the simulation data as observations. Using simulations, Welsch et al. (2007) have tested the robustness of the DAVE algorithm, the line-of-sight (LOS) magnetogram equivalent of the DAVE4VM algorithm (Schuck 2008), and found a good agreement between the simulated and the DAVE-derived velocities.

3.2.3. Current properties

Current-carrying magnetic fields are understood to be the building block for understanding the flares and CMEs drivers. The current properties computed here (see Table 3) allow us to quantify the energy stored by the magnetic field relative to its lowest energy state given its boundary distributions, i.e., potential. The vertical current density component J_z and the associated moments can be estimated through the classical Ampère’s law. The total, net, positive, and negative currents, respectively, noted as I_tot, I_net, and I_±, are estimated using various reckonings through the AR area. The heterogeneity and chirality current density components are also derived, according to Zhang (2001) (5)with B = Bb the first term of equation (5) refers to the current of chirality and the second relates to the heterogeneity, perpendicular by construction to the magnetic field. Thus, in the case of the heterogeneity current dominates over the chirality component, the AR magnetic field is far from the force-free field hypothesis. The associated four moments for each component are also evaluated, as well as the total and net heterogeneity and chirality currents, , . We also computed the net current originating from each polarity, , as described in Leka & Barnes (2003a) (6)as well as the vertical contribution to the current helicity density h_c along the vertical axis. Indeed, the current helicity density is defined as B·J but since only the vertical current component J_z can be deduced from the observations, the current helicity computed here is thus limited to the vertical contribution alone. From this partial estimation of the helicity density, we then derived the total and net partial helicity, and , characterizing the current helicity imbalance over the ARs. Bao et al. (1999) examined the relation between flare activity and the vertical contribution of the current helicity, and found that the time variation of partial current helicity is higher in flaring in ARs, suggesting partial and most probably total current helicity as a valuable eruptive indicator.

Determining whether or not the net electric current is neutralized over individual polarities of ARs is crucial for some theoretical flare and CMEs models (e.g., Melrose 1991; Parker 1996; Titov & Démoulin 1999; Forbes 2010). If the AR currents are fully neutralized, the net current integrated over one photospheric polarity is set to zero. Recent observations and simulations tend to confirm the existence of non-neutralized ARs (e.g., Török et al. 2014, Dalmasse et al. 2015, and references therein). Forbes (2010) argued that neutralized ARs may inhibit the eruption process. Hence, the measurement of direct and return current within ARs could provide insight into the flaring activity. The computations of both the direct current I_d and the return current I_r have been carried out following Török et al. (2014), who used the same simulations as this present work. They carried on the computation at the base of the corona, while we only concentrated on the surface, i.e., photospheric measurements. Several combinations of I_d and I_r have also been tested: I_d + I_r, |I_d/I_r| as well as the same quantities normalized by the total flux Φ_tot of the AR (see Sect. 3.2.1).

3.2.4. Lorentz forces

The Lorentz forces are thought to be an important parameter for solar eruptions. In particular, Fisher et al. (2012) argued that in eruptive flare ARs, the impulse coming from the Lorentz forces is dominant over all other forces, suggesting that an observed change in the Lorentz force before the eruption could be a pre-eruptive flare signature. To test this proposition in our simulations, we used the expressions derived by Fisher et al. (2012) in order to compute the horizontal and vertical components F_x, F_y, F_z of the Lorentz forces. We also estimated the total Lorentz force F and the normalized horizontal and vertical components , , (Table 4).

3.2.5. Magnetic polarity inversion line (MPIL) properties

Within an AR, the magnetic field is sheared, stressed, and twisted, thus deviating from the potential state of minimum energy. In particular, near the magnetic polarity inversion line (MPIL), the magnetic field can be almost perpendicular to the potential field (see e.g., Hagyard et al. 1990; Falconer et al. 2003). This region of strong magnetic shear is particularly well observed in X-rays, corresponding to the sigmoidal structures with an overall shape of “S” or an inverse “S” for bipolar ARs. Evaluating the properties of the magnetic field and currents near the MPIL provides a quantitative measure of the AR’s overall non-potentiality, even though some studies, such as Li et al. (2005), did not find obvious strongly sheared magnetic neutral line or strong overlying arcade associated with a X3 flare event. Since the degree of non-potentiality seems correlated with eruptive flare productivity (Falconer et al. 2002), the MPIL and near-MPIL areas’ properties are potentially proficient predictors.

We computed a series of parameters relying on MPIL lengths, most of them introduced by Falconer et al. (2003, 2006, 2008) (Table 5). The total length of the strong field neutral line L_s is defined as the length of the MPIL where is greater than a given threshold . The quantity L_ss defined the length of the MPIL for which the observed horizontal magnetic field B_h is greater than the given threshold and the shear angle is greater than . The total strong-gradient length of the MPIL L_sg is specified for MPIL regions where is greater than the previously defined threshold and where the horizontal gradient of the vertical magnetic field .We introduce here a new measure of the AR non-potentiality through the total strong-current length L_sc, measured along the MPILs, where and the current density J_z is greater than a given threshold .

The quantities WL_ss, WL_sg were first introduced by Falconer et al. (2008) while we define the new parameter WL_sc, measuring the current along the MPIL. WL_ss,WL_sg represent the integral of, respectively, the shear angle and the horizontal magnetic gradient , along the MPIL strong horizontal magnetic field portions. The WL_sc parameter corresponds to the current density J_z-integral along the MPIL portions of strong B_h. For the WL_ss and WL_sc parameters, only the MPIL portions where are taken into account, while WL_sg considers only the regions where . The two parameters L[Ψ > 80°] and L[ψ > 80°] from Leka & Barnes (2003a) are also tested, corresponding to the portions of the MPILs where, respectively, the 3D and the horizontal shear angles Ψ and ψ are greater than 80°. Finally, the unsigned flux near the MPIL is computed according to Schrijver (2007), given characteristic values around log R ~ 5.

Falconer et al. (2008) chose the following thresholds of G,  = 45° and G Mm⁻¹ while we initially adopt lower thresholds for G and G Mm⁻¹ given the small size and flux of our simulated ARs (see Sect. 3.1). For the L_sc parameter, we imposed a threshold  = 12 m Am⁻². It is worth noting that the choice of these thresholds can have substantial impact on the detection of eruptive flare signatures as we will show in Section 7, where a parametric study is proposed to estimate their influence.

4. Eruptive indicators evolution with B _mask = 30 G

For each output magnetogram of each simulation, we have computed the 99 parameters described in Section 3.2. For clarity, in this section, we only present the evolution of the most representative and interesting parameters among them. The parameters shown in Figures 3–5 have all been mentioned in previous studies as potential good indicators for eruptive and non-eruptive flaring activity. Even though the analysis method used in all these aforementioned works may differ from ours, i.e., using superposed epoch or forecasting methods, we used these studies as guidance to present our results and avoid quantitative comparisons. The displayed time window starts at t = 50 t ₀, once the flux tube has significantly emerged, the emergence starting at t = 35 t ₀. Below t = 50 t ₀, none of the parameters exhibits a significant variation, and for clarity, we thus restrained the displayed time window from t = 50 t ₀ to t = 150 t ₀, allowing us to analyze both prior and post eruptive flare (occurring at t ~ 120 t ₀) physical conditions. In the following, we do not attempt to provide physical interpretation of the full variation of the parameters, which is beyond the scope of this paper; rather we focus on the detection of significant changes in one or more parameters, according to the eruptive nature of the simulations only and occurring at a reasonable stage prior to the eruption.

Fig. 3. Some parameter evolution for the seven simulations, as a function of time in t 0 units (t 0 = 55.5 s). From (A–F), the following predictor evolution is respectively displayed: the total current helicity , the total current I tot, the mean value of the free energy , the total flux Φ tot, the total Lorentz force F and the Schrijver R value. The warm colors indicate the evolution for the three eruptive simulations {WD E, MD E, SD E}, respectively, corresponding to the magenta, orange, and red solid lines. Conversely, the evolution of the non-eruptive simulations {ND, WD, MD, SD} is displayed using the cold colors, respectively, as follows: purple, deep blue, light blue, and green solid lines. The eruption time is indicated by the vertical dashed gray line.

Figure 3 displays the evolution of some physical properties (see caption for details). The total partial current helicity (Fig. 3A), the total current I _tot (Fig. 3B), the free magnetic energy (Fig. 3C), and the total Lorentz force F (Fig. 3E) are in the top-ranked eruptive indicators by Bobra & Couvidat (2015), whereas the total unsigned flux ϕ _tot (Fig. 3B) and have been ranked in the best predictors for 48 hours flare forecasting by Bobra & Ilonidis (2016). The R value (Fig. 3F) was empirically established as a powerful eruptive indicator by Schrijver (2007). The parameters , I _tot, ϕ _tot, F, and the R value exhibit a very similar evolution, possibly due to an existing correlation between each parameter (Barnes et al. 2016). There is no evidence of physical changes prior to the eruptive flare, at about t ~ 120 t ₀ (vertical dashed gray line) for the eruptive simulations.

On the other hand, the mean free magnetic energy displays different behavior for the set of simulations, with a different peak value at t = 70 t ₀. The magnetic free energy is highly dependent on the region of interest selected using the B-mask described in Section 3.1. This selected area increases more rapidly than the free magnetic energy per surface unit after t = 70 t ₀, due to the emergence process. Therefore, the mean excess magnetic energy is observed to slowly decrease, but this is an artifact engendered by the size of the AR area. Using a constant mask as a test, the free magnetic energy is increasing as expected during the formation of a flux rope. However, the pattern does not correlate with the flare-eruptive nature of the simulations (warm vs. cold colors solid lines), but rather with the coronal arcade magnetic field strength. Higher peak values correspond to weaker coronal arcade, although that for a given coronal arcade magnetic field strength, the non-eruptive simulation shows systematically slightly higher values. Thus, the behavior of the indicator is not only dependent on the eruptive nature of the simulations, but is also coronal-field-strength dependent. This kind of behavior is not useful for flare forecasting, as we look for clear changes in some physical properties between stable and non-stable simulations only, in a way that some thresholds can, for example, be imposed. The property changes do not allow such an approach, since the same free magnetic energy peak could correspond to both eruptive and non-eruptive ARs, with different coronal arcade magnetic strengths.

Figure 4 presents the evolution of some magnetic field properties, such as geometry and current (see caption for details). The time evolution of the kurtosis of the twist parameter κ(α) is displayed in Figure 4A. The kurtosis of the horizontal magnetic field κ(B _h ) (Fig. 4C) is among the more efficient predictors in the four-variable discriminant analysis of Leka & Barnes (2003b), applied on flaring/non-flaring ARs. The mean gradient of the horizontal magnetic field (Fig. 4F) has been found by Bobra & Ilonidis (2016) to be the best predictive variable for 24 h flare forecasting, while the mean shear angle (Fig. 4D) is also in the top-ten best variables for both the 24 and 48 h predictions of the same study. We also examine the predictive capabilities of the relative magnetic helicity time variation (Fig. 4E) and the direct current I _d (Fig. 4B) since both have been suggested to play an important role in eruptive flare mechanisms (see Sect. 3.2 and e.g., Nindos & Andrews 2004; Dalmasse et al. 2015, and references therein).

Fig. 4. Same as Figure 3, but for different parameters. From (A–F), the following indicator evolutions are respectively displayed: the kurtosis of the twist parameter κ(α), the direct current , the kurtosis of the horizontal magnetic field κ(B h ), the mean value of the shear angle , the time variation of the relative magnetic helicity dH m/dt and the mean value of the vertical magnetic field B z .

As before, none of these series of parameters exhibits a clear eruptive signature. The direct current I _d and the relative helicity variation show very similar evolutions, even after the eruption at t ~ 120 t ₀. The kurtosis of the twist parameter κ(α) shows slightly higher values for the SD simulations, but apart from that difference, the evolution is almost the same for the whole simulation set. The κ(B _h ) presents a strong peak at t = 60 t ₀, with distinct values as a function of the simulations, but as for the , this behavior does not depend on the eruptive nature. Instead, the peak intensity depends on the strength of the coronal arcade, with stronger coronal magnetic field associated with stronger κ(B _h ) peak. The horizontal gradient also presents a strong peak, at about t ~ 66 t ₀, whose magnitude depends once again on the overlying field strength. As the coronal field gets smaller, the horizontal gradient increases, no matter the stable/unstable nature of the simulation. Notwithstanding, the mean value of the shear angle displays a slightly distinct behavior between eruptive and non-eruptive simulations. The parameter exhibits somewhat greater values for the eruptive numerical experiments, notably for the MD E and SD E simulations compared to that of the stable simulations. Both the eruptive and the non-eruptive remain stable once the flux tube has significantly emerged, i.e., after t = 70 t₀. However, there is no behavioral change before and after the eruption at t = 120 t₀, a required feature for a parameter to be a reliable predictor. As the eruption cannot be detected observing only the parameter, the mean shear angle is not a fully discriminant parameter.

Figure 5 shows the evolution for six of the MPIL properties extracted from the magnetogram series. All of them have been established as a powerful measure of the non-potentiality of ARs by Falconer et al. (2003, 2006, 2008), and are therefore potentially good predictors for both eruptive and non-eruptive flare activity. For this class of parameters, clear eruptive flare signatures are observable, and the evolution of eruptive and non-eruptive simulations is significantly different, except perhaps for the length of the strong-shear MPIL L_ss (Fig. 5A). For each of these parameters, a rise is undeniably observed during the early stage of the emergence for the three eruptive simulations, i.e., between t = 60 t₀ and t ~ 95 t₀. During the whole emergence process, the parameter evolutions of the SD E and MD E simulations have a significantly different comportment relative to the others, while those of the WD E simulation become similar to the non-eruptive simulations once the eruption is imminent, after t ~ 95 t₀. For the strong-shear length L_ss, both the SD E and MD E L_ss lengths are significantly longer, by respectively of factors 2 and 3, than that of the stable simulations, whereas the WD E L_ss length is barely 1.5 longer and only at the very early stage of the emergence, over a short time range. As a result, the L_ss is not a very strong predictor, since the eruptive signature quickly disappears, at least for weak coronal magnetic field. The parameter WL_sg (Fig. 5D) shows a similar trend, with a peak at t ~ 80 t₀ only about 1.2 times that of the non-eruptive simulations, making it a poorly robust predictive parameter.

Fig. 5. Same as Figure 3, but for different parameters. From (A–F), the following indicator evolutions are respectively displayed: the length of the strong-shear MPIL Lss, the integral of shear angle along the MPIL WLss, the length of the strong-gradient MPIL Lsg, the integral of the Bz horizonal gradient along the MPIL WLsg, the length of the strong-current MPIL Lsc and the integral of the current along the MPIL WLsc.

However, the other parameters, WL_ss (Fig. 5B), L_sg (Fig. 5C), L_sc (Fig. 5E), and WL_sc (Fig. 5F), are clearly more robust predictors, showing a distinct enhancement prior to the eruption, even for the WD E simulations. The WL_ss parameter is 3–8.5 times greater than that of the stable simulations, while the L_sg length is 2–3 times higher, but on a longer time range, between t = 75 t₀ and t = 95 t₀. The L_sc parameter presents a sharp peak at t = 80 t₀, between 2.8 and 5.4 times greater than the non-eruptive values. The WL_sc MPIL property is the most efficient eruptive indicator, since the three eruptive simulations reach similar values at the beginning of the emergence process, being about 8 times greater than the non-eruptive WL_sc measurements. For this latter parameter, the influence of the coronal arcade is reduced, at least during the early emergence stage.

All these MPIL properties provide comprehensible eruptive flare signatures in these models, prior to the eruption. Apart from the L_ss and WL_sg parameters, for which the signature exists but is weak, the L_sg, WL_ss, L_sc, and WL_sc parameters provide robust measurements of pre-flare eruptive conditions. As the MPIL properties are correlated with the magnetic configuration of the flux emergence, these parameters provide measurements of the magnetic complexity of the ARs. Since the eruptive flare formation and ejection are related to the quadrupolar versus bipolar nature of the emergence for this simulation set, a clear eruptive flare signature can be detected. The newly defined L_sc and WL_sc MPIL-current properties, as well as the WL_ss parameter appear to be the better predictor, with significant changes between eruptive and non-eruptive simulations. However, none of the other 93 parameters, related to current, magnetic field geometry, or properties, provides unambiguous eruptive flare signatures, since no significant changes have been detected in their evolution.

5. Influence of data masking

The eruptive indicators are highly sensitive to the AR area selected to perform the computation, as seen in Section 4, where the flux emergence process makes the AR area increase very rapidly. Bobra & Couvidat (2015) already highlighted that the AR parameters are highly sensitive to the data masking, and pointed out that a study to optimize such parameters is necessary. In this section, we present, as an example, the same results as in Section 4, but this time using a B_mask = 100 G. This new threshold leads to a reduced AR area, and because most of the parameters are AR-area dependent, most of the results are affected.

A variety of data masking methods can be found in previous studies. The automated system Solar Monitor Active Region Tracker (SMART) for detecting ARs and their associated properties using the Solar and Heliospheric Observatory (SOHO)/Michelson Doppler Imager (MDI; Scherrer et al. 1995) uses a static B_z threshold of 70 G to remove the background (Higgins et al. 2011). The global features provided by the Space-Weather HMI Active Region Patches (SHARPs; Bobra et al. 2014) pipeline select only the pixels for which the magnetic field strength is above the disambiguation threshold of about 150 G (see also Hoeksema et al. 2014). Using the ground-based data from the University of Hawai’i Vector Magnetograph (Mickey et al. 1996), Leka & Barnes (2003a) used only pixels above 3σ detections. Falconer et al. (2008) imposed the B_z component to be greater than 150 G, while Falconer et al. (2011) used either a threshold of 25 or 35 G, both studies using the MDI data. Al-Ghraibah et al. (2015) also used a 3σ binary mask to exclude noisy MDI data, while Mason & Hoeksema (2010) imposed a static threshold of 100 G.

Figure 6 displays the same parameter evolutions as Figure 5, but excluding pixels for which the total magnetic field B is lower than 100 G. The seven simulations exhibit the same evolution over time, reaching similar values close to the eruption starting time. In comparison with the previous results using a lower mask threshold, the rise observed for the six MPIL-related parameters of the eruptive simulations is lost, and the characteristic eruptivity signature is no longer discernible. This is due to the higher threshold masking imposed to the initial data, excluding initially from the dataset a fraction of the MPIL, and therefore reducing the computation domain of the MPIL properties.

Fig. 6. Same as Figure 5, but using a mask threshold Bmask = 100 G.

In order to visualize the data masking impact on the detection of the MPIL, Figure 7 displays the B_z magnetograms for the SD (Figs. 7A–7C) and the SD E (Figs. 7D–7F) simulations, at t = 100 t₀. The color scale is saturated between −400 and 400 G to highlight the two polarities. Since the orientation of the arcade is opposite in SD versus SD E simulations (see Sect. 2 for details), the magnetic configuration is purely bipolar for the non-eruptive simulations, whereas the eruptive magnetic topology is quadrupolar. This can be clearly observed in Figures 7A and 7D, where the MPIL is only located between the two polarities for the stable SD simulation (Fig. 7A), whereas the MPIL possesses additional portions at the polarity external edges for the SD E simulation. However, as the masking threshold increases, the external MPIL is reduced, and for the case where B_mask = 100 G, both MPILs, either in the eruptive or non-eruptive simulation, are almost exactly the same. Consequently, the deviation between stable and unstable simulations observed in the MPIL properties for B_mask = 30 G (see Fig. 5) is significantly reduced with B_mask = 100 G (see Fig. 6). The eruptive and non-eruptive simulations cannot be distinguished in this case, and thus the MPIL indicators become inefficient. This shows how much the initial mask threshold B_mask should be carefully chosen, since it may be crucial to be able to detect significant eruptive signatures and therefore make reliable flare forecasting.

Fig. 7. Magnetograms of the SD (A–C) and SD E (D–F) simulations, including the MPIL (red lines). The MPIL width has been dilated by a factor 4 in order to increase its visibility. The masking threshold increases from left to right, with respective values of 30, 50, and 100 G.

6. Impact of the noise on the detection of pre-eruptive signatures

Noise is an important issue for observational analysis. There are many different sources of error affecting the photospheric magnetic field measurements, e.g., photon noise, detector noises, spacecraft radial velocity inducing periodic systematic uncertainties (Hoeksema et al. 2014), or uncertainties associated with the inversion of the Stokes parameters. In the present work, these different sources are not individually treated, which is beyond the scope of this paper. Instead, we simply include a random Gaussian noise to our data in order to evaluate noise impact on the detection of pre-eruptive flare signatures.

We used a Monte-Carlo scheme, randomizing the magnetograms using Gaussian perturbations, with a standard deviation of 3.5 G. The error associated with HMI magnetic field data is typically about 10 G for the line-of-sight (LOS) component (V. Bommier, private communication) and can be generally one order of magnitude larger for the horizontal components (Hoeksema et al. 2014). However, our characteristic magnetic field values are about 3–4 times smaller than those that are typically observed (see Sect. 2.2). Therefore, the standard value of 10 G is not representative for our time series magnetograms. Accordingly, we adopted the characteristic noise standard deviation of 3.5 G, in order to remain consistent between our simulated observations and the relative noise scale. For this preliminary study, we assume the same error for the three components. Noise levels of 1 G and 5 G have also been tested with no sensitive differences with the results presented here. All the parameters investigated in this work are then derived from these time series of noisy magnetograms, following the same approach. First, for each simulation and for each time step, 50 different noisy magnetograms are computed by randomizing the initial one. The indicators are computed 50 times from these 50 noisy magnetograms and then averaged together. The error associated with each parameter for each time step is assumed to be the standard deviation corresponding to the 50 computed noisy indicator series, while the mean value provides the simulated magnetic field measurements.

Figure 8 displays the same parameters as Figure 5 but including the Monte-Carlo estimation of the measurement errors. The overall behavior remains unchanged, although some specific changes occur for the L_sc and the WL_sc (Figs. 8E and 8F) indicators. The L_ss and the WL_ss predictors for the MD E and SD E simulations are still significantly higher than that of the non-eruptive simulations. However, for these two parameters, in the WD E simulation, the pre-eruptive signature is lost. The sharp increase of the L_sg and WL_sg indicators, observed between t = 60 and t = 100 t₀, is still observable, making them robust to noise perturbation. However, for the two current-related predictors L_sc and WL_sc, the impact of noise is more important, and for instance the pre-eruptive signature provided by the WL_sc indicator is completely dominated by the noise. For the L_sc quantity, the trend is different, and a peak at t ~ 75 t₀ is still detectable for the eruptive simulations. However, we observe a very different behavior for the WD and ND simulations compared to the clean data (see Fig. 4). This is due to complex effect of the noise, which generates additional pseudo-MPIL in strong-current regions.

Fig. 8. Same as Figure 5, but including the noise perturbation.

To summarize, the presence of noise may affect the detection of the potential pre-eruptive signature, with impact depending on the parameter considered. The L_sg and WL_sg parameters are mostly not affected, and their pre-eruptive sharp increases are still detectable. L_ss and WL_ss and L_sc are slightly affected, showing weaker pre-eruption peaks, and moderate effects for the WD E simulations. However, WL_sc is strongly affected and the pre-eruptive signature is no longer observable.

7. Parametric study of MPIL properties

As described in Section 4, given our controlled-case study, only the MPIL features are able to provide a clear eruptive signature. Seven of the parameters characterizing the MPILs, namely L_s, L_ss, L_sg, L_sc, WL_ss, WL_sg, and WL_sc depend on the four different thresholds , , and (see Sect. 3.2.5 for details). These thresholds all refer to the portion of the MPILs taken into account in the computation of the various MPIL properties. As for the data masking process used to isolate the AR core area, the calculation of the eruptive indicators is highly sensitive to these values. In order to optimize the choice of these criterion and quantify their influence, we present in this section a parametric study of the MPIL properties. The threshold values explored for each parameter have been chosen based on typical values commonly used in observational data. Since our simulations correspond to small AR, we explore the parameter space using smaller and characteristic values used in previous studies.

To illustrate how much the choice of the thresholds potentially affects the detection of the MPILs, Figures 9–11 represent the MPIL portions (displayed as white line) detected as a function of the threshold , , and for the SD (top rows) and SD E (bottom rows) simulations, in a similar fashion as Figure 7. The length of the MPIL (white line) depends on these different parameters. Figure 9 displays B_h maps on the background, and the parameter varies from left to right, while Figure 10 represents maps as a function of the threshold. Figure 11 displays J_z maps and the threshold changes from 0 (left) to 12 m Am⁻² (right). In this section, the objective is to start to investigate the effects of each thresholding parameter on the MPIL parameters that have been demonstrated to be efficient eruptive predictors, under the condition that the initial mask thresholding was low enough (see Sect. 4). In this section, we do not initially mask the data (see Sect. 3.1, i.e., B_mask = 0 G) as we focus on threshold impact. Obviously, it is worth noting that the conjugate effects of high mask threshold (see Sects. 3.1 and 5) and individual parameter thresholdings are even stronger on the detection of the eruptive signature and therefore highly worsen the results.

Fig. 9. Maps of the horizontal magnetic field Bh, for the SD (A–C) and the SD E simulation (D–F) at t = 100 t0. The white line represents the portion of the MPIL for which . The threshold increases from left to right, with values of  = 0, 50, and 100 G. No initial data masking is applied here, in order to highlight only the impact of the threshold on the computation of Lss, Lsc, WLss, and WLsc.

Fig. 10. Maps of the quantity, using the same layout as Figure 9. The white line now represents the portion of the MPIL for which the horizontal gradient of the magnetic field is higher than . From left to right the threshold increases from 0–100 G Mm−1.

Fig. 11. Maps of the vertical current Jz, using the same layout as Figure 9. The white line now represents the portion of the MPIL for which the current density . From left to right, the threshold increases, with respective values 0, 3, and 12 mA m−2.

The threshold has a strong effect on the detection of the whole MPIL: as seen in Figure 9, the quadrupolar configuration, corresponding to outer MPIL on the edge of the polarities (see Fig. 9D), is no longer detected as threshold is greater than 50 G, therefore corresponding to similar MPIL (white line in Figs. 9C–9F). As such, any MPIL properties computed using such a threshold will not provide reliable eruptivity signature. The threshold is less restrictive, and a small fraction of the MPIL outer portions (white line in Fig. 9F) is still detected for the SD E simulations. However, it is worth noting that the MPIL represented in Figure 10 is dilated to increase the visibility; actually the MPIL is thinner and only a few additional pixels are detected. On the other hand, the thresholding does not affect that much the larger measurements of eruptive MPIL: as seen in Figure 11, the external portions of the MPIL are still observed, even using a threshold of 12 mA m⁻². From these series of maps as a preliminary coarse evaluation of thresholding effects on the MPIL properties, we can conclude that the current thresholding is the most flexible parameter, since the quadrupolar configuration of the eruptive simulations can still be detected even when increasing threshold.

To further investigate the influence of the physical thresholds on the MPIL predictors, we compute the variations of the L_ss, L_sg, L_sc, WL_ss, WL_sg, and WL_sc indicators as a function of their associated threshold. The three first quantities depend on two different thresholds, making therefore their evolution even more sensitive to the chosen cutting values. Figure 12 displays the L_ss variations for both eruptive SD E (solid lines) and non-eruptive SD (dashed lines) simulations, as a function of the and thresholds. The is increasing from panel to panel, while the is fixed for a given color line. If the is below 25 G (Figs. 12A and 12C), we can still observe longer L_ss lengths for the SD E simulation, whatever the threshold. In these three cases, the SD E L_ss parameter is respectively about 4, 3, and 2 times longer than that of the SD simulation, a parameter difference still measurable. However, for higher , this difference is reduced and to discriminate between eruptive and non-eruptive simulations becomes difficult. On the other hand, the influence of the threshold is weak and whatever the chosen value, SD E and SD simulations can still be distinguished.

Fig. 12. Parametric evolution of the total length of the strong-shear MPIL Lss, as a function of the two thresholds and . The varies for each figure, with respective values of 0, 10, 25, 50, 75, and 150 G, from (A–F). For each panel, the Lss quantity is represented for the strong arcade eruptive (SD E; solid lines) and strong arcade (SD; dashed lines) simulations. The threshold is changed for each curve, using the values of 30° (blue lines), 45° (red lines), and 60° (magenta lines). The eruption time is denoted by the vertical gray dot dashed line.

Figure 13 is the same as Figure 12, but exploring the variations of the L_sg predictor as a function of and thresholds. As before, each panel corresponds to a given threshold, and each color line corresponds to a given . For this parameter, the influence of the threshold is lower, while the impact of the is significantly stronger. For example, for the given G Mm⁻¹ (red lines on each panel), whatever the imposed, the SD E L_sg length is still longer than that of the SD simulation. However, as the increases, the divergence between the SD E and SD curves rapidly disappears, making the eruptive signature undetectable. For Figures 12A and 12B, it is worth noting that the length L_sg (dashed blue line) of the SD simulation remains constant over the whole numerical experiment. This is due to the low imposed, below the photospheric arcade magnetic field mean value, conjugated to the G Mm⁻¹ threshold, allowing therefore to measure the whole MPIL on the entire numerical domain. Since SD is a non-eruptive simulation, the magnetic field configuration is bipolar, and the MPIL is almost kept constant.

Fig. 13. Same as Figure 12 but for the parameter Lsg, as a function of the two thresholds and . The varies in the same way as for Figure 12. The threshold is changed for each curve, using the values of 0 (blue lines), 10 (red lines), 25 (magenta lines), 50 (green lines), 75 (yellow lines), and 100 G Mm−1 (black lines).

Figure 14 shows in the same way as Figures 13 and 12 the variations of the newly introduced L_sc parameter, as a function of its two associated thresholds and . For this MPIL property, the influence of both thresholds is important. As for the L_ss parameter, detecting the eruptive signature, i.e., a longer L_sc for the eruptive simulations, becomes harder as exceeds 25 G. For G, a longer L_sc for the SD E simulation is however still detectable, but the difference between eruptive and non-eruptive simulations is less obvious. For higher (Figs. 14E and 14F), discriminating the simulations is impossible, whatever the current . On the other hand, increasing the also has a strong effect on the detection of the eruptive signature, even though small differences between eruptive and non-eruptive simulations persist. For instance, if G (Figs. 14B) and  = 3 mA m⁻² (blue lines), L_sc is almost four times longer before the eruption for the SD E simulation than that of the SD. If we increase to 24 mA m⁻² (yellow lines), the L_sc SD E to SD ratio decreases to 1.7.

Fig. 14. Same as Figure 12 but for the parameter Lsc, as a function of the two thresholds and . The threshold is changed for each curve, using the values of 3 (blue lines), 6 (red lines), 12 (magenta lines), 18 (green lines), and 24 mA m−2 (yellow lines).

From these three plot analyses, general dependence trends can be deduced. The L_ss parameter is strongly dependent on the threshold, and for greater than 50 G, the eruptivity signature is no longer measurable. However, the influence of the shear-angle threshold is rather weak, and whatever the threshold assumed, L_ss parameters still allow us to discriminate between the SD E and SD simulations. Conversely, the L_sg threshold is weakly dependent on the threshold, while the gradient threshold has a strong impact on the detection of eruptive ARs. For G Mm⁻¹, the distinction between both AR types is hardly observable. Finally, the L_sc quantity is impacted by the choice of both its associated thresholds, even if small remnants of eruptive signature are still detectable for high thresholding values.

The WL_ss, WL_sg, and WL_sc parameters only depend on the threshold (see Table 5). Figure 15 displays the variations of these indicators for both the SD E (solid lines) and the SD (dashed lines) simulations, varying the threshold, in order to optimize the choice of this parameter. Figure 15A displays the WL_ss parameter evolutions, and as the threshold increases, the WL_ss parameter decreases for the SD E simulation while it remains stable for that of the SD non-eruptive. Using a low excluding value allows us to detect the external portion of the MPIL (blue, red, and magenta lines) and therefore to discriminate between the eruptive and non-eruptive simulations, while thresholds larger than 50 G (green, yellow, and black lines) remove this eruptive signature, characterized by a larger WL_ss for eruptive ARs.

Fig. 15. Parametric evolution of the WLsc (A), WLsg (B), and WLss (C) parameters, as a function of the threshold . Solid lines represent the parameter evolutions for the SD E simulation, while dashed lines are computed for the SD stable simulation. The varies from 0 (blue lines) to 150 G (black lines), with the following intermediate values 10 (red lines), 25 (magenta lines), 50 (green lines), and 75 G (yellow lines).

Figures 15A and 15B present respectively the WL_sc and WL_sg and parametric analysis. For these two parameters, the thresholding impact is weak, and even using high threshold allows to distinguish between eruptive and non-eruptive simulations. The WL_sg property (see Fig. 5D) is about 2–2.5 times higher for the SD E simulation than that of the SD, whatever the adopted threshold imposed on . The newly defined WL_sc parameter shows sharp variations for the eruptive simulation, increasing quickly by a factor about 4, whereas the non-eruptive evolution remains smooth over the whole emergence process. These sharp evolutions are observed whatever the B_h threshold imposed on the horizontally observed magnetic field, demonstrating that the noise-free WL_sc is a robust predictor for detecting flare producing ARs.

From these parametric studies of the MPIL properties WL_ss, WL_sg, and WL_sc, clear conclusions can be drawn. The WL_ss predictor is highly dependent on the threshold process, and the eruptive signature, corresponding to higher WL_ss in eruptive ARs, is rapidly lost as increases. Consequently, if the WL_ss is able to detect imminent eruption in an AR under certain conditions (i.e., low masking threshold B_mask, see Sect. 5), this may not be the most reliable eruptive indicator to be used for flare forecasting, unless different thresholds are simultaneously tested. On the other hand, both noise-free WL_sg and WL_sc parameters appear as very robust eruptive predictors since the choice of the threshold does not affect the detection of the eruptive signature.

8. Summary and conclusions

The predictability of magnetic properties has been barely estimated using numerical simulations (see Sect. 1). In this work, we have studied the reliability of some eruptive flare indicators using the 3D parametric MHD simulations of flux emergence in Leake et al. (2013, 2014), leading to the self-consistent formation of either stable or unstable flux ropes. A set of four stable and three eruptive simulations have been used, corresponding to the partial emergence of a twisted magnetic flux tube into a stratified atmosphere, where a coronal overlying field is present. The eruption is triggered by a combination of external magnetic reconnection between the dipole and the overlying magnetic field, followed by internal reconnection allowing the vertical expansion and thus the ejection of the flux rope. The behavior of the emergence has been investigated through a range of initial coronal arcade strengths and orientations. Depending on the coronal arcade orientation, the magnetic configuration is either dipolar or quadrupolar, leading to respectively stable or unstable coronal flux rope (see Sect. 2.1).

Parameters have been examined in order to detect whether some physical changes specific to eruptive flaring ARs could be detected and measured, thus establishing a certain relation to flare occurrence. We rigorously analyzed the simulations as observations, including data masking and magnetic flux transport velocity derivation from the DAVE4VM code (Schuck 2008), using time series of 2D-plane magnetograms extracted from the 3D numerical datasets. A list of 99 parameters has been tested, including indicators currently used for both eruptive and non-eruptive flare operational forecasting as well as new quantities, such as helicity or current-weighted neutral line lengths L_sc and WL_sc.

From the 99 predictors’ list used in this work, only the parameters relying on MPIL properties (see Table 5) demonstrated predicting eruptive flare capabilities. The other parameters showed a very similar evolution for both simulations, and no eruptive signature was detectable (see Sect. 4). The strong-shear MPIL length L_ss, the strong-gradient MPIL length L_sg and the strong-current MPIL length L_sc are significantly longer for the unstable simulations due to the quadrupolar magnetic configuration generating additional external MPIL (see e.g., Fig. 7). The length increase rate depends on the coronal magnetic field strength: stronger coronal dipole is associated to longer property-weighted MPIL lengths. However, a strong rise of these three lengths is systematically observed during the early stage of the emergence, whatever the arcade field strength. The WL_ss, WL_sg, and the WL_sc parameters present very similar trends, making them also promising eruptivity predictors.

The detection of eruptive flare signatures, i.e., higher MPIL properties for eruptive ARs, is dependent on the initial masking process. Indeed, to isolate the ARs from the background magnetic field, a thresholding mask is usually applied when observations are analyzed. Notwithstanding, we demonstrated that the choice of the initial masking threshold is crucial for detecting MPIL property variations. The same analysis has been presented for both B_mask = 30 and 100 G (see Sect. 5), showing that the eruptive flare signatures disappear using the higher mask threshold. Noise in the measurements can also scramble the eruptivity signatures: through a Monte-Carlo scheme, we estimated the noise influence by including random perturbation to our time series magnetograms. Our results show that the four L_ss, L_sg, WL_ss, and WL_sg parameters are not strongly impacted by noise and their associated peaks prior to the eruption are still detectable. However, the L_sc and WL_sc are more strongly impacted: the increase of L_sc is weaker, although still measurable, but the WL_sc parameter no longer allows to discriminate between eruptive and stable simulations.

In addition, the MPIL properties, associated with valuable eruptive flare predictabilities, depend not only on the initial masking process and the noise, but also on additional thresholds imposed on physical properties, such as current, magnetic field, or gradients. We also investigated the impact of these physical thresholds through a parametric study measuring the detectability of eruptive signature as a function of the thresholding process (see Sect. 7). Results show that the three MPIL lengths L_ss, L_sg, and L_sc and the WL_ss parameter are strongly sensitive to the choice of the physical threshold, whereas WL_sg and WL_sc are robust relative to threshold changes.

Our study is limited in terms of AR size and complexity that we cannot explore given the computational cost of such simulations. Given the small scale of our MHD simulations, we do not catch the same distribution of magnetic field strengths as observations. Because of the small flux (in the 10²¹ Mx range), and the absence of interaction between granulation and magnetic field, the intermediate strength fields are not caught. Hence, by applying the B-mask, which is done at these intermediate strengths, some features from the simulations may be lost. Besides, eruptive flare indicators have only been tested for a given AR flux and size, corresponding to the smallest observed flaring AR class (see Sect. 2.2 for details). Still, this type of analysis is yet relevant to connect flare physical models and observations, and provide a comprehensive perspective of what can be done using actual observations. Future work will be extended to AR simulations with different sizes and fluxes.

We therefore conclude that from our 99 predictors’ list, WL_sg and L_sg are the best eruptive flare indicators for these model tests. The other parameters relying on MPIL properties tested in this study, namely L_ss and L_sc, WL_sc and WL_ss should be tested using various physical thresholds. The current-related parameters seem to be more sensitive to noise, even if a more detailed analysis of uncertainty sources is needed. Apart from the physical thresholds, we also recommend the testing of various masking processes for actual observational analysis, using both low and high values in order to detect complex MPILs, potentially indicating an imminent eruptive flare activity.

Acknowledgments

This work has received partial support by the European Unions Horizon 2020 Research and Innovation Programme under Grant Agreement No. 640216 (Flare Likelihood and Region Eruption Forecasting [FLARECAST]). JEL is supported by NASA’s Heliophysics program. Simulations used in this study were carried out using a grant of computational time from the DoD’s High Performance Computing Program. The editor and the authors thank two anonymous referees for their assistance in evaluating this paper and improving its quality.

References

All Tables

All Figures

Fig. 1. Emergence of the convection magnetic flux tube into the solar atmosphere at t = 110 t 0 for both the stable SD (A) and the unstable SD E (B) simulations. The cyan lines trace the coronal arcade magnetic field, originating from the lower boundaries, whereas the red lines belong to the magnetic flux tube and originate from the y = ±maxy side boundaries. The gray scale indicates the magnetic field strength at the surface.

In the text

Fig. 2. Magnetograms of the Bz component at t = 100 t0 for both the MD (A) and MD E (B) simulations. Both non-masking (left) and masking (right) magnetograms are displayed, in order to show the reduction of the region of interest for further analysis. Eruptive and non-eruptive simulations present very similar magnetograms, except near the external edge of the polarities, due to the quadrupolar geometry above the surface.

In the text

Fig. 3. Some parameter evolution for the seven simulations, as a function of time in t 0 units (t 0 = 55.5 s). From (A–F), the following predictor evolution is respectively displayed: the total current helicity , the total current I tot, the mean value of the free energy , the total flux Φ tot, the total Lorentz force F and the Schrijver R value. The warm colors indicate the evolution for the three eruptive simulations {WD E, MD E, SD E}, respectively, corresponding to the magenta, orange, and red solid lines. Conversely, the evolution of the non-eruptive simulations {ND, WD, MD, SD} is displayed using the cold colors, respectively, as follows: purple, deep blue, light blue, and green solid lines. The eruption time is indicated by the vertical dashed gray line.

In the text

Fig. 4. Same as Figure 3, but for different parameters. From (A–F), the following indicator evolutions are respectively displayed: the kurtosis of the twist parameter κ(α), the direct current , the kurtosis of the horizontal magnetic field κ(B h ), the mean value of the shear angle , the time variation of the relative magnetic helicity dH m/dt and the mean value of the vertical magnetic field B z .

In the text

Fig. 5. Same as Figure 3, but for different parameters. From (A–F), the following indicator evolutions are respectively displayed: the length of the strong-shear MPIL Lss, the integral of shear angle along the MPIL WLss, the length of the strong-gradient MPIL Lsg, the integral of the Bz horizonal gradient along the MPIL WLsg, the length of the strong-current MPIL Lsc and the integral of the current along the MPIL WLsc.

In the text

	Fig. 6. Same as Figure 5, but using a mask threshold Bmask = 100 G.
In the text

	Fig. 7. Magnetograms of the SD (A–C) and SD E (D–F) simulations, including the MPIL (red lines). The MPIL width has been dilated by a factor 4 in order to increase its visibility. The masking threshold increases from left to right, with respective values of 30, 50, and 100 G.
In the text

	Fig. 8. Same as Figure 5, but including the noise perturbation.
In the text

	Fig. 9. Maps of the horizontal magnetic field Bh, for the SD (A–C) and the SD E simulation (D–F) at t = 100 t0. The white line represents the portion of the MPIL for which . The threshold increases from left to right, with values of  = 0, 50, and 100 G. No initial data masking is applied here, in order to highlight only the impact of the threshold on the computation of Lss, Lsc, WLss, and WLsc.
In the text

	Fig. 10. Maps of the quantity, using the same layout as Figure 9. The white line now represents the portion of the MPIL for which the horizontal gradient of the magnetic field is higher than . From left to right the threshold increases from 0–100 G Mm−1.
In the text

	Fig. 11. Maps of the vertical current Jz, using the same layout as Figure 9. The white line now represents the portion of the MPIL for which the current density . From left to right, the threshold increases, with respective values 0, 3, and 12 mA m−2.
In the text

Fig. 12. Parametric evolution of the total length of the strong-shear MPIL Lss, as a function of the two thresholds and . The varies for each figure, with respective values of 0, 10, 25, 50, 75, and 150 G, from (A–F). For each panel, the Lss quantity is represented for the strong arcade eruptive (SD E; solid lines) and strong arcade (SD; dashed lines) simulations. The threshold is changed for each curve, using the values of 30° (blue lines), 45° (red lines), and 60° (magenta lines). The eruption time is denoted by the vertical gray dot dashed line.

In the text

	Fig. 13. Same as Figure 12 but for the parameter Lsg, as a function of the two thresholds and . The varies in the same way as for Figure 12. The threshold is changed for each curve, using the values of 0 (blue lines), 10 (red lines), 25 (magenta lines), 50 (green lines), 75 (yellow lines), and 100 G Mm−1 (black lines).
In the text

	Fig. 14. Same as Figure 12 but for the parameter Lsc, as a function of the two thresholds and . The threshold is changed for each curve, using the values of 3 (blue lines), 6 (red lines), 12 (magenta lines), 18 (green lines), and 24 mA m−2 (yellow lines).
In the text

Fig. 15. Parametric evolution of the WLsc (A), WLsg (B), and WLss (C) parameters, as a function of the threshold . Solid lines represent the parameter evolutions for the SD E simulation, while dashed lines are computed for the SD stable simulation. The varies from 0 (blue lines) to 150 G (black lines), with the following intermediate values 10 (red lines), 25 (magenta lines), 50 (green lines), and 75 G (yellow lines).

In the text