Solar Magnetic Flux Rope Eruption Simulated by a Data-Driven Magnetohydrodynamic Model

TL;DR

This paper presents a data-driven magnetohydrodynamic model that simulates the eruption of solar magnetic flux ropes, aligning well with multi-wavelength observations and advancing understanding of solar eruption mechanisms.

Contribution

It introduces a novel data-driven MHD simulation approach using observed magnetic and velocity fields to accurately model flux rope eruptions.

Findings

The model reproduces magnetic flux rope evolution during eruptions.

Simulation results match observed flare ribbons and flux rope morphology.

A fixed boundary condition yields similar eruption simulations.

Abstract

The combination of magnetohydrodynamic (MHD) simulation and multi-wavelength observations is an effective way to study mechanisms of magnetic flux rope eruption. We develop a data-driven MHD model using the zero-$\beta$ approximation. The initial condition is provided by nonlinear force-free field derived by the magneto-frictional method based on vector magnetic field observed by the Helioseismic and Magnetic Imager (HMI) aboard the Solar Dynamics Observatory (SDO). The bottom boundary uses observed time series of the vector magnetic field and the vector velocity derived by the Differential Affine Velocity Estimator for Vector Magnetograms (DAVE4VM). We apply the data-driven model to active region 11123 observed from 06:00 UT on 2011 November 11 to about 2 hours later. The evolution of the magnetic field topology coincides with the flare ribbons observed in the 304 and 1600 \AA $ $ wavebands by the Atmospheric Imaging Assembly. The morphology, propagation path, and propagation range of the flux rope are comparable with the observations in 304 \AA . We also find that a data-constrained boundary condition, where the bottom boundary is fixed to the initial values, reproduces a similar simulation result. This model can reproduce the evolution of a magnetic flux rope in its dynamic eruptive phase.

Figures

Equations

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{v})={\nu }_{\rho }{\nabla }^2(\rho -{\rho }_{\mathrm{i}}),$$

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla \cdot (\mathbf{v}\rho \mathbf{v}-\mathbf{BB})+\nabla (\frac{{\mathbf{B}}^2}{2})=\mu \nabla \cdot [2\mathcal{S}-\frac{2}{3}(\nabla \cdot \mathbf{v})\mathcal{I}],$$

$$\frac{\partial \mathbf{B}}{\partial t}+\nabla \cdot (\mathbf{vB}-\mathbf{Bv})=-\nabla \times (\eta \mathbf{J}),$$

$$\mathbf{J}=\nabla \times \mathbf{B},$$

$$T=\{\begin{array}{ll}{T}_0& 0.0\le h<h_0\\ k_T(h-h_0)+T_0& h_0\le h<h_1\\ T_1& h_1\le h<13.8\end{array},$$

$$\frac{\mathrm{d}p}{\mathrm{d}h}=-g\rho ,$$

Full text

Sun: activity — Sun: corona — Sun: coronal mass ejections (CMEs) — Sun: flares — Sun: magnetic fields — Sun: photosphere

1 Introduction

Solar activity, such as prominence/filament eruptions, coronal mass ejections (CMEs) and flares are mostly related to magnetic flux rope eruptions, where the flux rope either pre-exists before or is formed during eruption (Ouyang et al., 2017). These eruptions expel magnetized plasma, with a huge amount of energy and helicity, and high energy particles into the interplanetary space and possibly disturb the space environment around the Earth. It has been proposed that one important driving mechanism for magnetic flux rope eruptions is due to ideal MHD instability, such as the torus instability (equivalent to the loss of equilibrium) for isolated flux ropes (Démoulin & Priest, 1988; Forbes & Isenberg, 1991; Lin & Forbes, 2000; Kliem & Török, 2006; Démoulin & Aulanier, 2010; Kliem et al., 2014), or tilt-kink instability for parallel current systems (Keppens et al., 2014), or coalescence-kink instability in anti-parallel current systems (Makwana et al., 2018). At the same time, the resistive process of magnetic reconnection could form the magnetic flux rope at the buildup phase, or facilitates the eruption in the driving phase (Antiochos et al., 1999; Moore et al., 2001; Aulanier et al., 2010; Chen, 2011). However, the interplay and feedback between ideal and resistive processes are still not fully understood. In this sense, it is crucial to build a working model for studying and, hopefully, predicting magnetic flux rope eruptions.

The buildup phase of a magnetic flux rope in a solar active region can be modeled with a sequence of magnetic equilibria, and has been widely studied by the nonlinear force-free field (NLFFF) models (e.g. Canou et al., 2009; Savcheva & van Ballegooijen, 2009; Guo et al., 2010; Inoue et al., 2012; Jiang et al., 2014; Yang et al., 2016). The dynamic process of magnetic flux rope eruptions can only be studied by magnetohydrodynamic (MHD) simulations. Here, we focus on a special class of such simulations, namely, data-driven and data-constrained MHD simulations, which are defined as simulations initiated and/or driven by observed data. Various data-driven numerical simulations differ mainly in three aspects. In the aspect of initial conditions, most of the data-driven and data-constrained MHD models are implemented with the initial condition provided by NLFFF models, derived by either the magneto-frictional (Kliem et al., 2013; Jiang et al., 2016; Inoue et al., 2018) or the Grad-Rubin method (Amari et al., 2018). In the aspect of the boundary conditions, they could be provided from a time series of observations, which is called data-driven simulations (Jiang et al., 2016), or be assigned numerically with fixed values or zero-gradient extrapolations, which is called data-constrained simulations (Kliem et al., 2013; Inoue et al., 2018; Amari et al., 2018). In the aspect of the MHD model, either the zero-$\beta$ model (Kliem et al., 2013; Inoue et al., 2018; Amari et al., 2018), which omits the gas pressure, gravity, and the energy equation, or the ideal adiabatic MHD model (Jiang et al., 2016) has been used thus far.

The data-driven and data-constrained MHD simulations have shown their strength in studying solar eruptions. For instance, we could study the buildup, triggering, and driving processes of a magnetic flux rope and its eruption; or investigate the roles of ideal and resistive processes in a flux rope eruption; or predict a flux rope eruption. Our motivation in this paper is to develop a data-driven MHD model and to directly confront it with observations. The observations of the vector magnetic field and extreme-ultraviolet (EUV) images are introduced in Section 2. The MHD model is described in Section 3. The results are presented in Section 4. Summary and discussions are provided in Section 5.

2 Observations

We use the vector magnetic field, UV, and EUV observations in active region NOAA 11123 on 2010 November 11 as the input and benchmark for the data-driven MHD simulations. The vector magnetic field is obtained from the spectro-polarimetric observations provided by the Helioseismic and Magnetic Imager (HMI; Scherrer et al., 2012; Schou et al., 2012; Hoeksema et al., 2014) aboard the Solar Dynamics Observatory (SDO; Pesnell et al., 2012). Its cadence is 12 minutes and spatial sampling is $0.5^{\prime\prime}$ per pixel. The line-of-sight component of the vector magnetic field in active region 11123 at 06:00 UT is shown in Figure 1a. Active region 11123 is a newly emerging region trailing behind the long-living preceding active region 11121. Mandrini et al. (2014) studied in detail the magnetic topology and eruptive activities in this active region complex. Galsgaard et al. (2015) studied the mechanisms of up-flows in this active region using a data-driven MHD simulation. In this study, we focus on the C4.7 class flare starting at 07:16 UT and peaking at 07:25 UT on 2010 November 11 accompanied with a small CME with a speed of about 250 km s*-1* (Schmieder et al., 2014). The Atmospheric Imaging Assembly (AIA; Lemen et al., 2012) aboard SDO provides continuous observations in three UV-visible wavebands and seven EUV wavebands with a cadence of 12 seconds and spatial sampling of $0.6^{\prime\prime}$ per pixel. We use the UV and EUV images, especially in the 1600, 304, and 171 Å wavebands, observed by SDO/AIA to study the eruptive process and as the benchmark to test the MHD models. The pre-eruptive coronal loops and filament are shown in the 171 Å and 304 Å wavebands at 06:00 UT in Figure 1c and 1d, respectively.

3 MHD Model

The lower corona is in a low-$\beta$ condition, where $\beta$ is the ratio between the gas pressure to the magnetic pressure. Thus, we adopt a zero-$\beta$ MHD model to simulate the dynamics in the corona. This model omits the gradient of gas pressure and gravity, and it also omits the energy equation. Only three physical variables, namely, the density, velocity, and the magnetic field, need to be solved:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $\rho$ is the density, $\mathbf{v}$ is the velocity, $\mathbf{B}$ is the magnetic field, $\nu_{\rho}$ is the density diffusion coefficient, $\rho_{\mathrm{i}}$ is the density at $t=0$, $\mu$ is the dynamic viscosity coefficient, $\mathcal{S}_{ij}=\frac{1}{2}(\partial v_{i}/\partial x_{j}+\partial v_{j}/\partial x_{i})$ is the strain rate tensor, $\mathcal{I}$ is the unit tensor, and $\eta$ is the resistivity. Equations (1)–(4) are in the dimensionless form such that the vacuum permeability is assumed to be $\mu_{0}=1$. In the numerical setup, all the physical quantities and parameters are normalized by their corresponding typical factors, namely, $L_{0}=1.0\times 10^{9}$ cm, $t_{0}=85.9$ s, $\rho_{0}=2.3\times 10^{-15}\ \mathrm{g\ cm^{-3}}$, $\mathbf{v}_{0}$ = $1.2\times 10^{7}\ \mathrm{cm\ s^{-1}}$, and $\mathbf{B}_{0}=2.0$ G. The MHD equations are solved by the open source Message Passing Interface Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al., 2003, 2012; Porth et al., 2014; Xia et al., 2018).

The source terms on the right hand side of Equations (1)–(3) are included with different purposes. Since the gradient of the gas pressure is omitted, the density evolution cannot be correctly handled without a density diffusion term in Equation (1). This term is used to smooth the density distribution and maintain the initial density stratification in the evolution and the density diffusion coefficient is finally set to be $\nu_{\rho}=0.008$ by trial and error. The viscous term in Equation (2) is adopted for numerical stability and the dynamic viscosity coefficient is set to be $\mu=0.05$. The resistive term in Equation (3) is used to control the explicit magnetic diffusion and study the effect of magnetic reconnection. We use $\eta=0$ in the first two ideal cases and non-zero values in other cases.

The initial condition for the density $\rho$ is provided by a stratified atmosphere profile to simulate the solar atmosphere from the photosphere to the corona. We assume a stepwise function for the temperature:

[TABLE]

where $T_{0}=0.006$, $T_{1}=1.0$, $h_{0}=0.35$, $h_{1}=1.0$, and $k_{T}=(T_{1}-T_{0})/(h_{1}-h_{0})$ in dimensionless units, or equivalently, $T_{0}=6.0\times 10^{3}$ K, $T_{1}=1.0\times 10^{6}$ K, $h_{0}=3.5\times 10^{8}$ cm, and $h_{1}=1.0\times 10^{9}$ cm in cgs units, which closely mimics the actual temperature-density variation with height. The density profile is derived by solving the hydrostatic equation:

[TABLE]

where the gas pressure $p=\rho T$ in the dimensionless form. The density, $\rho$, on the bottom is chosen to be $1.0\times 10^{8}$ (namely, $2.3\times 10^{-7}$ g cm*-3*), and it drops to 1.1 at $h=13.78$, i.e., the top of the computation box. Although we prescribe the density in this realistic fashion, the actual zero-$\beta$ simulation does not have pressure or temperature information, which is why we use the density diffusion in Equation (1). The initial condition for the velocity is zero for all the three components, i.e., $\mathbf{v}=0$ at $t=0$.

The initial condition for the magnetic field $\mathbf{B}$ is provided by the NLFFF model derived by the magneto-frictional method (Guo et al., 2016a, b). We use the vector magnetic field by SDO/HMI at 06:00 UT on 2010 November 11 as the boundary condition for constructing the NLFFF model. The series name for the vector magnetic field in the Joint Science Operations Center is “hmi.B_720s”. Two more processing steps in addition to the SDO/HMI pipeline have been applied to the vector magnetic field. One is the correction of projection effects using the formula in Gary & Hagyard (1990). The other is preprocessing (Wiegelmann et al., 2006) to remove the net Lorentz force and torque in order to conform to the force-free assumption. The projected and preprocessed vector magnetic field is shown in Figure 1b. We use a vertically stretched grid to resolve the steep gradient of the density profile. The grid cell sizes from the bottom to top are geometric series. Therefore, the cell size in the vertical direction is $2.77\times 10^{-3}$ on the bottom and 0.37 on the top, while the cell sizes on the horizontal directions are uniformly $0.76\times 10^{-1}$. The computation box in the range $[x_{\mathrm{min}},x_{\mathrm{max}}]\times[y_{\mathrm{min}},y_{\mathrm{max}}]\times[z_{\mathrm{min}},z_{\mathrm{max}}]=[-20.47,-6.79]\times[-35.57,-21.89]\times[0.1,13.78]$ is resolved by $180^{3}$ cells. The location and field of view of the computation box are shown in Figure 1a. The finally derived NLFFF is displayed in Figure 1c and 1d, which shows the comparisons between magnetic field lines with flare loops and a filament. We note that the computation box has been projected to the heliocentric coordinate system to compare with observations directly. The formulae to do the back projection are described in Guo et al. (2017).

We need to specify the boundary conditions for density, velocity, and magnetic field for each of the six boundaries ($x_{\mathrm{min}},x_{\mathrm{max}},y_{\mathrm{min}},y_{\mathrm{max}},z_{\mathrm{min}},z_{\mathrm{max}}$). There are two ghost layers on each boundary to conform with the requirements of the numerical scheme, which is three-step time integration, HLL Riemann solver, and Koren limiter space reconstruction. To test the effects of the boundary conditions, we prepare two different cases:

Case I

data-driven boundary condition. The density is fixed to be the initial value on all the six boundaries. The velocity is fixed to be zero on the top and four side boundaries, and is set to be the data derived by the Differential Affine Velocity Estimator for Vector Magnetograms (DAVE4VM; Schuck, 2008) on the inner ghost layer of the bottom boundary, and zero-gradient extrapolation on the outer ghost layer of the bottom boundary. Since the velocity can only be derived at discrete times with a cadence of 12 minutes, the needed data with much higher cadence for the bottom boundary condition are computed with a linear interpolation in time. Finally, the magnetic field is provided by zero-gradient extrapolation on the top and four side boundaries, and is set to be the observed data on the inner ghost layer of the bottom boundary, and zero-gradient extrapolation on the outer ghost layer. Similar to the velocity, linear interpolation is used to fill the data gap between each magnetic field observation with a cadence of 12 minutes.

We take some measures to guarantee the boundary condition conforming with the induction equation and $\nabla\cdot\mathbf{B}=0$ condition. First, we add a diffusive term $\delta\nabla(\nabla\cdot\mathbf{B})$ in the induction equation (Keppens et al., 2003) to guarantee the solenoidal condition in the computation domain. Second, the normal component of the vector magnetic field on the outer ghost layer at the bottom boundary (and the two ghost layers at the other five boundaries) is reset by requiring $\nabla\cdot\mathbf{B}=0$ to guarantee the solenoidal condition for the boundaries. Third, the velocity field is derived by DAVE4VM, a method that is fully consistent with the induction equation (Schuck, 2008).

Case II

data-constrained boundary condition. The velocity on the bottom boundary (both ghost layers) is fixed to be zero. The magnetic field on the inner ghost layer for the bottom boundary is fixed to be the initial observed data, namely, the vector magnetic field observed at 06:00 UT on 2010 November 11. The others are the same as Case I.

4 Results

The evolution of the magnetic field in the process of the flux rope eruption is shown in Figure 2. At 06:00 UT (Figure 2a), a magnetic flux rope as indicated by the sheared and twisted magnetic field lines lies under three magnetic null points (Mandrini et al., 2014), one of which is shown by the surrounding spine-fan shaped field lines. We also compute the quasi-separatrix layers (QSLs) to highlight the boundaries of different magnetic domains using a recently developed three-dimensional (3D) QSL computation method (Pariat & Démoulin, 2012; Yang et al., 2015). Figure 2b reveals that the flux rope is under the fan QSL. As time goes on, the two ends of the flux rope curve up (Figure 2c), and it stretches and distorts the fan QSL (Figure 2d). The flux rope is finally detached from the bottom and erupts into the corona as shown in Figure 2e in a sample snapshot at 07:04 UT. The QSLs surrounding and interleaving the flux rope also break the fan QSL (Figure 2f). At about 07:26 UT, the front of the flux rope rises to a high altitude and reaches the side boundary of the computation box (Figure 2g). The QSLs associated with it develop a very complex structure distributed over the boundary and body of the flux rope (Figure 2h). The propagation path of the front of the flux rope is not along the radial direction, but is inclined to the southeast.

Figure 3 shows the comparison of the QSLs on the bottom with the flare ribbons observed by SDO/AIA in both 304 and 1600 Å wavebands. There are three major QSLs on the bottom as shown in Figure 3a at 07:22 UT, labeled the eastern QSL E, the middle QSL M, and the western QSL W, respectively. There are also three flare ribbons appearing in the 304 Å image at 07:22 UT. It is found that Ribbons E, M, and W coincide spatially with parts of QSLs E, M, and W, respectively (Figure 3b). Two flare ribbons (E and M) appear in 1600 Å, which also coincide with parts of the corresponding QSLs (Figure 3c). Flare ribbons always appear at the footpoints of QSLs. But QSLs could exist even in potential magnetic field with no electric current and free magnetic energy. Only when particles are accelerated or plasma are heated in QSLs, would flare ribbons appear at the footpoints of these QSLs. We note that when the QSLs are compared with the 304 and 1600 Å observations, they have been back projected to heliocentric coordinates.

The evolution of the flux rope simulated by the zero-$\beta$ MHD model is also compared with SDO/AIA 304 Å observations. The results in Figure 4 indicate that the MHD simulations conform with the observations in three aspects, namely, in morphology, eruption path, and propagation range. First, Figure 4a shows a snapshot of the 304 Å observation at 07:22 UT, where the erupting bright strand as indicated by the white arrow is generally believed to be a flux rope. The flux rope is highly curved, whose height is apparently larger than the distance between its two footpoints. Some flare loops under the erupting flux rope are so bright that their intensity is saturated in the 304 Å image. The curved flux rope and the underlying cusp-shaped flare loops can also be found in the MHD model as shown in Figure 4b. Second, the inclined eruption path is reproduced by the MHD model. It is not along the solar radial direction but inclined to the southeast. The projection of the eruption path coincides with the 304 Å observations as shown in Figure 4c–4f and the animation attached to Figure 4. Finally, the flux rope in the observations erupted into interplanetary space and formed a CME (Schmieder et al., 2014). The animation attached to Figure 4 also indicates that the flux rope continues to erupt close to the boundary of the computation box.

Next, we compare quantitatively the eruption speeds of the flux rope both from the observation and from the MHD simulation. The speed observed in 304 Å is measured by the time-distance diagram as shown in Figure 5a, where the slice is indicated in Figure 4a. We pinpoint the front of the flux rope at each observation time. The measurement is repeated ten times to estimate errors. We also measure the 3D positions of the flux rope front at selected snapshots in the MHD simulation (Figure 5b). The projection of the 3D positions to the SDO viewpoint is shown as the green dots in Figure 5c. Similar to the time-distance profile in 304 Å, the profile of the MHD simulation also shows a slow rising stage followed by a rapid rising stage. The only difference is that the flux rope starts to rise about 20 minutes earlier in the MHD simulation than that in the 304 Å observation. Additionally, the flux rope rises a little bit faster and higher in the simulation than in the observation in the slow rising stage, while it is slower (30.1 km s*-1*) in the simulation than in the observation (83.8 km s*-1*) in the rapid rising stage. If we shift the time-distance profile in the simulation 20 minutes, the evolution processes are similar to each other as shown in Figure 5c.

Finally, we compute the decay index along the eruption path of the flux rope to provide a qualitative explanation of the time-distance profile. Since the flux rope does not propagate along the radial direction of the Sun, the actual eruption direction is determined by a nonlinear least-square regression of the 3D points of the flux rope fronts, whose unit vector is defined as $\mathbf{e}_{\mathrm{pro}}$. The decay index is computed by $n=-\mathrm{d}\log B_{\mathrm{ex,pol}}/\mathrm{d}\log{s}$, where $B_{\mathrm{ex,pol}}$ is the field strength of the poloidal component of the external magnetic field, $s$ is the distance along the propagation path. Here, we use the potential field at the initial time, 06:00 UT, to represent $B_{\mathrm{ex}}$, whose poloidal component is defined as follows. We require that the unit vector of the poloidal direction, $\mathbf{e}_{\mathrm{pol}}$, is perpendicular to both $\mathbf{e}_{\mathrm{pro}}$ and the polarity inversion line, which is approximately along $\mathbf{e}_{x}$ in this case. These conditions require that $\mathbf{e}_{\mathrm{pol}}\times(\mathbf{e}_{\mathrm{pro}}\times\mathbf{e}_{x})=0$. The field strength of the poloidal component of the external magnetic field is computed by the projection of the potential field to the poloidal direction, $B_{\mathrm{ex,pol}}(s)=\mathbf{B}_{\mathrm{pot}}(s)\cdot\mathbf{e}_{\mathrm{pol}}$. The decay index at $s$ is thus obtained. It is not displayed against the distance $s$ along the propagation path, but is the projected distance in the plane-of-sky as shown in Figure 5d to compare with the observations directly.

It is found that the flux rope starts to rise rapidly at about 14.0 Mm as shown by the green dots in Figure 5c. And at about 13.7 Mm (Figure 5d) the decay index crosses the canonical critical value of 1.5. We note that the critical decay index of torus instability is not strictly 1.5 but ranges from 1 to 2 depending on different current paths (straight, semicircular or others), which has been found by either theoretical or numerical methods (Bateman, 1978; van Tend & Kuperus, 1978; Kliem & Török, 2006; Török & Kliem, 2007; Démoulin & Aulanier, 2010; Fan, 2010; Olmedo & Zhang, 2010). Taking the range of 1 to 2 and referring to Figure 5d, it is found that the torus unstable distance ranges from about 13.0 to 14.0 Mm, which is a small range and is consistent with the turning position (14.0 Mm) where the eruption velocity changes from slow to fast. This result is in agreement with the prediction of the torus instability. From Figure 5c, we find that the flux rope starts to rise at about 5.0 Mm as shown by the green dots. This distance is far below the torus unstable region as found in Figure 5d using the decay index range of 1 to 2. Therefore, the flux rope needs additional physical mechanisms to rise to the torus unstable region. Magnetic reconnection in null points or below the flux rope (similar to the tether-cutting reconnection) could perform such a role to raise a flux rope to the torus unstable region.

5 Summary and Discussion

We have developed a data-driven MHD model using the zero-$\beta$ MHD equations, which are solved by the open-source code MPI-AMRVAC (Keppens et al., 2003, 2012; Porth et al., 2014; Xia et al., 2018). This model is applied to active region 11123 that was observed by SDO/HMI and SDO/AIA. The initial condition of the MHD simulation is provided by the NLFFF model constructed from the vector magnetic field at 06:00 UT on 2010 November 11 and the magneto-frictional method (Guo et al., 2016a, b). The boundary condition on the bottom is based on vector magnetic field observations and vector velocity derived by the DAVE4VM method (Schuck, 2008). The data-driven MHD simulation reproduces the eruption process of the magnetic flux rope. The QSLs in the vertical direction show that the flux rope lies under a complex spine-fan structure, which contains three magnetic null points (Mandrini et al., 2014). Parts of the QSLs on the bottom surface coincide with the flare ribbons observed in 304 and 1600 Å. The morphology, propagation path, and propagation range of the flux rope eruption in the MHD simulation resemble those in the SDO/AIA 304 Å observations. The height-time profiles of the flux rope front in both the MHD simulation and the 304 Å observations have a two stage evolution, one of which is slow and the other is rapid. These similarities demonstrate physics-based numerical predictions of solar eruptions. However, it is noted that the flux rope in the MHD simulation erupts about 20 minutes earlier than the observation, and the velocity of the flux rope is slower in the simulation than that in the observation. These discrepancies indicate the limitations of the present zero-$\beta$ simulation, which asks for further improvements, such as including the energy equation and considering the effects of thermal conduction and radiative cooling.

To test the effects of different boundary conditions, we prepare a data-constrained boundary, namely Case II as described in Section 3. The simulation results of Case II are very similar to Case I in terms of morphology, propagation path and range of the flux rope. This similarity reveals two noteworthy aspects of the MHD model. First, magnetic flux rope eruptions possess varying temporal and spatial scales. Compared to the dynamic eruptive phase, the buildup phase of a magnetic flux rope is relatively long, which may last for several days or even weeks. The data-driven effects should be more important in this long phase than in the eruptive one. Therefore, the data-driven boundary condition is indispensable in modeling the buildup phase, while it might be negligible close to an eruption. Second, this result implies that we could predict the eruption of a magnetic flux rope in principle. Since the effect of the boundary condition could be neglected when initial conditions are close enough to the dynamic eruptive phase, we might only use data in advance of the eruption to predict this following evolution. Furthermore, the apparent indifference of our simulation between the data-driven and data-constrained approaches confirms once more that ideal MHD instabilities, such as the torus instability in this particular configuration, unavoidably lead to eruptions. The initial rise to build up the torus unstable state in the data-constrained simulation might be caused by the unbalanced Lorentz force in the numerical NLFFF model, or a finite but small resistivity in the numerical diffusion of the simulation. Since it is well-known that subtle details in internal pitch or external field variations play a decisive role in determining the growth rate of ideal MHD instabilities (see, e.g., Goedbloed & Poedts, 2004), it remains challenging to get precise agreements in full eruption dynamics.

We note that the above two cases are only influenced by numerical resistivity since the explicit resistivity was set to be 0 so far. We also investigate the effects of tuning the resistivity, $\eta$, in Equation (3). Two sets of tests are performed. In the first set, $\eta$ is uniform in the whole computation box. It is found that when $\eta\leq 5\times 10^{-5}$, the simulation result is very similar to that with $\eta=0$. When $\eta\geq 2\times 10^{-4}$, we find noticeable differences in the simulation result, where the rising speed of the erupting flux rope is slower and the regions with large electric current start to change drastically. The observations do not show such changes. The experiments imply that the resistivity, which is derived from our inherent numerical diffusion in the previous $\eta=0$ runs, is of an acceptable level to obtain CME eruptions that are comparable with observations. In the second set, $\eta=\eta_{0}$ when $J<J_{c}$, and $\eta=\eta_{2}[(J-J_{c})/J_{c}]^{2}+\eta_{0}$ when $J\geq J_{c}$, where $\eta_{0}$ and $\eta_{2}$ are the resistivity of the background and the coefficient of anomalous resistivity, respectively, and $J_{c}$ is the critical electric current density. In a reasonable range of these parameters, the results do not show better evolution of the magnetic field compared to observations. A possible reason is that the excitation of anomalous resistivity should not be determined by the current density only. To simulate the effect of magnetic reconnection in accordance with observations, we need to specify its location and timing with more constraints or resort to more advanced coupled MHD and particle in cell (PIC) methods, as recently demonstrated on coalescing islands in Makwana et al. (2018).

The authors thank the anonymous referee for constructive comments. The SDO data are available by courtesy of NASA/SDO and the AIA and HMI science teams. YG thanks Chaowei Jiang for useful discussions on the boundary condition. YG, MDD, and PFC are supported by NSFC (11773016, 11733003, and 11533005), NKBRSF 2014CB744203, the fundamental research funds for the central universities 020114380028, and Jiangsu 333 Project (BRA2017359). CX is supported by NSFC (11803031). RK is supported by KU Leuven (GOA/2015-014). The numerical calculations in this paper have been done using the computing facilities of the High Performance Computing Center (HPCC) in Nanjing University.