MHD simulation of Solar Eruption from Active Region 11429 Driven by Photospheric Velocity Field

As driven by solar eruptions, the solar-terrestrial environment often experiences variations, which are known as space weather, and forecasting the space weather precisely is not only an important scientific topic but can also avoid damage of the sensitive on-ground and space-based critical infrastructures. Though many theoretical models have been proposed and significant process has been made in understanding the triggering mechanism of solar eruptions (Forbes, 2000; Chen, 2011; Schmieder et al., 2013; Priest, 2014), reproducing the whole life-span from quasi-static stage to eruption using numerical models constrained and driven by observed vector magnetogram possess unprecedented capabilities in revealing the onset mechanism of the real eruption events, and can potentially be used for accurate space weather forecast (Jiang et al., 2022; Jiang, 2022).

Previous study reproduced the whole process of energy accumulation and release successfully (Jiang et al., 2016), showing an MHD system can be driven to erupt by inputting time series of vector magnetograms at the bottom boundary (B-driven). There are also other data-driven models, in which the evolution of MHD system is driven by the electric field (E-driven, e.g., Cheung & DeRosa 2012; Hayashi et al. 2018; Pomoell et al. 2019; Price et al. 2019) or the velocity field (V-driven, e.g., Hayashi et al. 2019; Guo et al. 2019; Liu et al. 2019; He et al. 2020; Zhong et al. 2021) on the photosphere (bottom boundary). Though the B-driven method can fully match the magnetogram, it will introduce considerable errors of magnetic divergence from the bottom boundary. This shortage will vanish in E-driven model, however, deriving both the induction and potential components of the electric field on the photosphere is not an easy task (Fisher et al., 2010, 2012, 2015, 2020), and the photospheric flow also needs to be properly set to follow the Ohm’s law. In most of the theoretical models of solar eruption, the key structure in favor of eruption is assumed to be formed through the movement of the footpoints of magnetic field lines (Moore & Labonte, 1980; Moore & Roumeliotis, 1992; Antiochos et al., 1999; Lin & Forbes, 2000; Jiang et al., 2021b), which is driven by the horizontal flow, and the vertical component of the photospheric flow will be responsible for the flux emergence process. Therefore, with the photospheric velocity field determined, the photospheric magnetic field can be generated self-consistently. Furthermore, in the V-driven approach, there is no need to solve the complex momentum equation at the bottom boundary (which is the most time-consuming part in solving MHD equations). Due to these advantages, the V-driven method has attracted many previous studies to focus on this topic. For example, with the velocity field derived by DAVE4VM method (Schuck, 2008), Hayashi et al. (2019) used the projected normal characteristics method to update the physical variables other than velocity at the bottom boundary. However, the total magnetic energy kept almost the same level of that of the initial state without obvious magnetic energy injection. Jiang et al. (2021a) updated the magnetic field by solving directly the magnetic induction equation at the bottom boundary and their model can inject the magnetic energy from bottom boundary successfully. He et al. (2020) drove the magnetic evolution by inputting the DAVE4VM velocity field and vector magnetogram simultaneously. The formation process of an magnetic flux rope (MFR) was obtained in their model. Unfortunately, these simulations didn’t drive the system to erupt. The only work we know that obtained an eruption by the velocity field derived from observation was shown in Kaneko et al. (2021). The velocity field was derived from the electric field on the photosphere using Ohm’s law and then was input at the bottom boundary in their zero-beta model to drive the system to erupt. Two eruptions they produced were identified from the evolution curves of the magnetic and kinetic energies, however, the kinetic energy showed an overall increase without obvious quasi-static evolution during the pre-eruption stage. Before the first eruption, the kinetic energy was comparable with its peak during the first eruption and the amount of the release of the magnetic energy was also too small, which didn’t show the feature of a typical eruption event, i.e., a large amount of magnetic energy was converted into kinetic energy impulsively. As described above, the whole energy accumulation process from quasi-static evolution (of typically tens of hours) to impulsive eruption has not been realized in a self-consistent way using the V-driven MHD model, and this is one of the motivations of this work.

In this Letter, we applied a V-driven model to investigate the evolution and eruption of a well-studied active region (AR) NOAA 11429. Previous studies found persistent shearing flow and flux cancellation near the main polarity inversion line (PIL) of this AR (Shimizu et al., 2014; Zheng et al., 2017), which were suggested to be responsible for the eruptions on 2012 March 7, 9 and 10 (Dhakal et al., 2018, 2020). An analysis of the MFR reconstructed from vector magnetograms using the nonlinear force-free field (NLFFF) model suggests that the homologous eruptions are triggered by torus instability (TI) of the MFR (Chintzoglou et al., 2015). Zhang et al. (2021) suggested the helical kink instability may also take effect. Nevertheless, whether these mechanisms were at work requires to be further studied using dynamic modeling of this AR evolution and eruption, which is absent in all the previous studies and is the other motivation of this work. In our simulation, the dynamic process from the beginning of 2012 March 4 to the eruption on March 5 in this AR was reproduced self-consistently as driven by the photospheric velocity derived from vector magnetogram using DAVE4VM method. Our simulation shows that an MFR was generated near the main PIL before the flare onset and the initiation of this eruption event depended mainly on TI of the preformed MFR.

AR 11429 showed a complex $\beta\gamma\delta$ configuration and is very flare-productive, which has produced 3 X-class flares from 2012 March 5 to 7. It first appeared on the eastern solar limb on March 4 and was located on the eastern part of the solar disk before March 8. During this period, the AR kept developing as characterized by the increasing total unsigned magnetic flux, indicating the obvious flux emergence by vertical flow on the photosphere. Since it was the first X-class flare of this AR, here we focus on the initiation process of the X1.1 flare (which is accompanied with a halo CME moving at a speed of 1531 $\rm km~{}s^{-1}$) around 04:00 UT on March 5 as shown in the white box in Figure 1A. Before the flare onset, there was persistent shearing flow near the main PIL (Figure 1B) and as a result, the horizontal magnetic field there was highly sheared (Figure 1C), which indicates that a large amount of free magnetic energy is stored ready for an eruption. An hot loop first erupted away as shown in AIA 94 $\rm\AA$ (Lemen et al., 2012) at around 03:31 UT (Figure 1D and E) and after that a pair of hook shape flare ribbons appeared near the main PIL in AIA 1600 $\rm\AA$ at 03:36 UT (Figure 1F), i.e., the flare event started.

To understand the formation of the pre-eruptive coronal magnetic field and the triggering mechanism of this eruption, we used the DARE-MHD model (Jiang et al., 2016) to study the dynamic evolution of this AR. For saving the computational time, the strength of magnetic field from the magnetograms were reduced by a factor of 25 before being input into our code. The initial plasma density was set as a hydrostatic isothermal model with a fixed temperature as $T=1\times 10^{6}$ K and a modified solar gravity (Jiang et al., 2021b) to get a plasma background that mimics the real environment in the solar corona basing on two key parameters, the plasma $\beta$ and the $\rm Alfv\acute{e}n$ speed (in particular, the minimum $\beta$ is $6.3\times 10^{-4}$ and maximum $\rm Alfv\acute{e}n$ speed $V_{\rm A}\sim 4800$ km s^-1 in the final equilibrium we obtained below). We chose to use the magnetograms of HMI SHARP data set (Schou et al., 2012; Bobra et al., 2014) at 00:00UT on March 4 to reconstruct the initial magnetic field since there were no obvious MFR at that time. This magnetogram was smoothed first using Gaussian smoothing with FWHM of 6 pixels and a NLFFF model from this smoothed magnetogram was extrapolated by our CESE-NLFFF-MHD code (Jiang & Feng, 2012, 2013). Since the code (like many other NLFFF codes) does not give a perfect force-free solution but with residual Lorentz forces, we input the extrapolated field into the DARE-MHD model, along with the initial background plasma, to let the MHD system relax until the kinetic and magnetic energies were almost unchanged, i.e., an MHD equilibrium was obtained, and the initial state was ready.

At the bottom boundary, we solve the magnetic induction equation to update the magnetic field with the velocity field derived by DAVE4VM method to update the magnetic field on the photosphere. The DAVE4VM velocity was strengthened by a factor of 13.7 (determined by the ratio of the time series magnetograms’ original time cadence as 720 seconds to the time cadence in our simulation as $0.5\times 105$ seconds, i.e., $\frac{720}{0.5\times 105}$) to speed up our simulation and thus the time scale of quasi-static evolution prior to the eruption onset is shorten by the same times. At the side and top boundaries, all the variables are extrapolated from the neighboring inner points with zero gradient along the normal direction of the boundary surface and the normal component of magnetic field is further modified by the divergence-free condition. The Powell-source terms and diffusion control terms was used to deal with the divergence error of the magnetic field as described in Jiang et al. (2010). We set the computational domain sufficiently large as $[-368,368]$ Mm in both $x$ and $y$ direction and $[0,736]$ Mm in $z$ direction with grid resolution varies from $1^{\prime\prime}$ to $8^{\prime\prime}$ using adaptive mesh refinement. The highest resolution is used mainly for the regions with strong magnetic gradients and current density, in particular, the current sheets (CSs). Explicit value of magnetic resistivity was not used in our simulation and the magnetic reconnection was controlled by the resistivity of the numerical method only, which mimicked the low-resistivity plasma better. As the total unsigned flux of this AR kept increasing before the flare onset, we energized the MHD system by full 3D DAVE4VM velocity field ($v^{D}_{x},v^{D}_{y},v^{D}_{z}$) (will be referred to as V3D simulation) and horizontal component ($v^{D}_{x},v^{D}_{y},0$) (V2D) respectively, and a comparison of the results for these two simulations will show the importance of flux emergence through the vertical flow on the photosphere.

We have carried out a full 3D MHD simulation of an X-class flare eruption event on 2012 March 5 in AR 11429 using the V-driven DARE-MHD model. An MHD equilibrium was obtained by relaxing the NLFFF reconstructed by the CESE-MHD-NLFFF code and was set as the initial condition. Then the initial state was driven to evolve by the DAVE4VM velocity field on the bottom boundary of our simulation box. The analysis of the quasi-static evolution stage before the eruption onset shows the gradual formation of an MFR above the main PIL. When the MFR first appeared, it was relatively low and then it grew up into the torus unstable region, after which TI was triggered. Then the energy conversion process was accomplished by reconnection in the CS that is formed due to the stretching effect of the erupting MFR. The images of SDO observation and the important physical quantities computed from the magnetograms are reasonably consistent with our result in the pre-flare magnetic structure, the morphology of flare ribbons, evolution curves of the magnetic and kinetic energies as well as the total unsigned magnetic flux. The time duration of the quasi-static evolution process before the simulated eruption is also very close to the actual time scale before the flare onset.

Nevertheless, our simulation does not reproduce accurately the evolution of magnetic field as shown in the observed magnetograms. As can be seen in the last panel of Figure 3A, part of the magnetic flux was transported to and concentrated at the edge of the AR, while this pileup was not observed in HMI magnetograms. One likely reason for this flux pileup is that, the DAVE4VM method only solves the normal component of the magnetic induction equation in the least-square sense (Lumme et al., 2019), which may not reproduce the velocity precisely at every point and can only be used to approximate the overall distribution and evolution of magnetic flux in the AR. In addition, the magnetic flux on the photosphere should be dispersed and dissipated by granular and supergranular convection (Wang et al., 1989) as well as small scale turbulent diffusion in practical situation. Therefore, without dealing properly with these effects, the flux pileup will be more obvious than observations and as a result the magnetic energy injection rate of ‘V3D’ run will be overall higher than the ‘OB’ one as shown in Figure 2A (comparing the black dashed line and the black solid line). This unrealistic flux pileup may also contribute to the overshoot of the magnetic energy increase (Figure 2A, the mismatch of the black dashed line and the green line before the eruption), because it results in a very large magnetic gradient at the bottom boundary, thus even the finest spatial resolution in our simulation was inadequate to capture the too high gradient during the time duration close to the eruption onset (i.e., $t\in[97,122]$ minutes). The insufficient grid resolution led to the considerable numerical error there, which thus makes the magnetic energy increase higher than the magnetic energy injection in ‘V3D’ simulation. The proper settings of numerical diffusion and grid resolution, along with the improved methods for deriving the photospheric flow will need to be considered in future works for a more accurate data-driven simulation of solar eruptions.

To summarize, our simulation shows that, besides inputting the time series of vector magnetogram, the numerical model of solar corona can be driven to erupt by inputting the time series velocity field at the bottom boundary, in which the driver, i.e., the bottom velocity field is derived from time series of vector magnetograms. The numerical model we established here shows the possibility of driving the evolution of solar corona using different physical variables, which offers a straightforward way for revealing the eruption mechanisms in real events. Thus, such model has a great potentiality in forecasting the onset time as well as the strength of solar eruptions and evaluating the quantitative impact on space weather.  

This work is jointly supported by the National Natural Science Foundation of China (NSFC 41731067, 42030204, 42174200), the Fundamental Research Funds for the Central Universities (grant No. HIT.OCEF.2021033), and Shenzhen Science and Technology Program (grant Nos. RCJC20210609104422048 and JCYJ20190806142609035). The computational work was carried out on TianHe-1(A), National Supercomputer Center in Tianjin, China, and we thank Jun Chen for his informative and helpful discussions.