Magnetic fields dominate the dynamics in the Sun's upper atmosphere, the solar corona. On the solar surface, i.e., the photosphere, magnetic fields are seen to change continuously; magnetic flux emergence brings new flux from the solar interior into the atmosphere, and meanwhile the flux is advected and dispersed by surface motions such as granulation, differential rotation, and meridional circulation. Consequently, the coronal field evolves in response to (or driven by) the changing of the photospheric field, and thus complex dynamics occur ubiquitously in the corona, including the interaction of newly emerging field with the pre-existing one, twisting and shearing of the magnetic arcade fields, magnetic reconnection, and magnetic explosions which are manifested as flares and coronal mass ejections (CMEs).

As we are still not able to measure directly the three-dimensional (3D) coronal magnetic fields, numerical modeling has long been employed to reconstruct or simulate the coronal magnetic fields based on different assumption of the magnetohydrodynamic (MHD) equations, such as the most-frequently used force-free field model [1], which has been developed for over four decades. However, the force-free field assumption is only valid for the equilibrium of the corona, and it cannot be used to follow a continual and dynamic evolution of the magnetic fields. In recent years, with vector magnetogram data in the photosphere measured routinely with high resolution and cadence, data-driven modeling is becoming a viable tool to study coronal magnetic field evolution, which can self-consistently describe both the quasi-static and the dynamic evolution phases [e.g., 2–10]. Still, due to the limited constraint from observation, data-driven models are developed with very different settings from each other [11]. For simplicity, some used the magneto-frictional model [4, 10], in which the Lorentz force is balanced by a fictional plasma friction force. As such, the magnetic field evolves mainly in a quasi-static way [12], although in some case an eruption can be reproduced but it evolves in a much slower rate than the realistic one [10], because the dynamic is strongly reduced by the frictional force. Some used the so-called zero-β model [7, 9, 13], in which the gas pressure and gravity are neglected. The zero-β model might fail when there is fast reconnection in the field, in which the thermal pressure could play an important role in the dynamics in weak-field region of magnetic field dissipation. Therefore, it is more realistic and accurate to solve the full MHD equations to deal with the non-linear interaction of magnetic fields with plasma.

To drive the full MHD model, one needs to specify all the eight variables (namely plasma density, temperature, and three components of velocity and magnetic field, respectively) in a self-consistent way at the lower boundary. Previously, with only the magnetic field obtained from observations, we employed the projected-characteristic method [2, 14–17] to specify the other variables according to the information of characteristics based on the wave-decomposition principle of the full MHD system [e.g., 5, 18]. Specifically, the full MHD equations are a hyperbolic system that can be eigen-decomposed into a set of characteristic wave equations (i.e., compatibility equations), which are independent with each other; on the boundary, these waves may propagate inward or outward of the computational domain because the wave speeds (the eigenvalues) of both signs generally exist at the boundary. If the wave goes out of the computational domain through the boundary, it carries information from the inner grid to the boundary, and thus the corresponding compatibility equation should be used to constrain the variables on the boundary surface. Thus, if there are five waves going out, and with three components of magnetic field specified by observed data, the eight variables are fully determined. Ideally and in principle, the projected-characteristic method is the most self-consistent one with inputting data that is partially available (for example, in our case, only the data of magnetic field). However, such conditions, i.e., exactly five waves going out, are often not satisfied in the whole boundary and other assumptions are still necessary. Even worse, in areas near the magnetic polarity inversion line (PIL) where the normal magnetic field component is small, the Alfvén wave information goes mainly in the transverse direction rather than the normal one, and the projected-characteristic approach may fail. Another shortcoming of the method lies in its difficulty in code implementation, which needs to perform eigen-decomposition and solve a linear system of the compatibility equations on every grid point on the boundary to recover the primitive variables from the characteristic ones.

In this paper we test another way of specifying the bottom boundary conditions, which uses the surface velocity derived by the DAVE4VM method [19] to drive the MHD model. The differential affine velocity estimator (DAVE) was first developed for estimating velocities from line-of-sight magnetograms, and was then modified to directly incorporate horizontal magnetic fields to produce a differential affine velocity estimator for vector magnetograms (thus called DAVE4VM). It is generally accepted that the coronal magnetic fields are line tied in the dense photosphere and are advected passively by the surface motions of the photosphere, such as the shearing, rotational, and converging flows [20]. Thus, with the surface velocity in hand, we can solve the magnetic induction equation on the bottom boundary to update self-consistently the magnetic field to implement such a line-tying boundary condition, which is a much simpler way than using the projected-characteristic method. To illustrate the approach, we take the solar AR 12158 as an example to simulate its two-day evolution from 2014 September 8 to 10 during its passage of the solar disk. This AR is of interest since it produced an X1.6 eruptive flare accompanied with a fast CME with speed of ~1300 km s⁻¹, which is well-documented in the literature [21–23]. The AR is well-isolated from neighboring ARs, thus is suited for our modeling focused on a single AR. Within a few days prior to this major eruption, the AR developed from a weakly sheared magnetic arcade into a distinct sigmoidal configuration, indicating a continual injection of non-potential magnetic energy into the corona through the photospheric surface motion. Indeed, prior to the eruption, the major sunspot of the AR showed a significant rotation of over 200° in 5 days [21]. Based on the sigmoidal hot coronal loop seen immediately before the flare, many authors have interpreted it as a pre-existing magnetic flux rope which erupted and resulted in the flare and CME [e.g., 21–23], while some non-linear force-free field extrapolations appear not to support this [24, 25]. Thus, a fully data-driven MHD simulation can provide valuable insight in addressing this issue by following the dynamic evolution of the coronal magnetic field, although the main objective of this paper is to describe the methods.

In the following we first describe our model equation in section 2. The data-driven simulation consists of three steps, and first we constructed an MHD equilibrium based on a single vector magnetogram observed for the start time of our simulation, which is described in section 3. Then in section 4 we calculated the surface flow field using the DAVE4VM code with the time series of vector magnetograms. Finally, we input the flow field in the model to drive the evolution of the MHD system, as described in section 5. Summary and discussion are given in section 6.

We numerically solve the full MHD equations in a 3D Cartesian geometry by an advanced conservation element and solution element (CESE) method [26]. Before describing the model equations in the code, it is necessary to specify the quantities used for non-dimensionalization. Here we use typical values at the base of the corona for non-dimensionalization as given in Table 1. In the rest of the paper all the variables and quantities are written in non-dimensionalized form if not mentioned specially.

In non-dimensionalized form, the full set of MHD equations are given as

where J=∇×B, ν is the kinetic viscosity, and γ is the adiabatic index.

Note that we artificially add a source term −ν_ρ(ρ − ρ₀) to the continuity equation, where ρ₀ is the density at the initial time t = 0 (or some prescribed form), and ν_ρ is a prescribed coefficient. This term is used to avoid a ever-decreasing of the density in the strong magnetic field region, which we often encounter in the very low-β simulation. It can maintain the maximum Alfvén speed in a reasonable level, which may otherwise increase and make the iteration time step smaller and smaller and the long-term simulation unmanageable. Specifically, this source term is a Newton relaxation of the density to its initial value by a time scale of

where τ_A = 1/v_A is the Alfvén time with length of 1 (the length unit) and the Alfvén speed vA=B/ρ. Thus, it is sufficiently large to avoid its influence on the fast dynamics of Alfvénic time scales. As a result, we used ν_ρ = 0.05v_A in all the simulation in this paper.

No explicit resistivity is used in the magnetic induction equation, but magnetic reconnection is still allowed through numerical diffusion when a current layer is sufficiently narrow such that its width is close to the grid resolution. In the energy (or temperature) equation, we set γ = 1 for simplicity, such that the energy equation describes an isothermal process¹. The kinetic viscosity ν will be given with different values when needed, which is described in the following sections.

We first constructed an initial MHD equilibrium based on a single vector magnetogram taken for time of 00:00 UT on 2014 September 8 by SDO/HMI. Such an equilibrium is assumed to exist when the corona is not in the eruptive stage, and is crucial for starting our subsequent data-driven evolution. The vector magnetogram is preprocessed to reduce the Lorentz force and is further smoothed to filter out the small-scale structures that are not sufficiently resolved in our simulation. Reduction of the Lorentz force is helpful for reaching a more force-free equilibrium state [e.g., 27], and smoothing is also necessary to mimic the magnetic field at the coronal base rather directly the photosphere, because the lower boundary of our model is placed at the base of the corona [28]. The preprocessing is done using a code developed by Jiang et al. [29], which is originally designed for NLFFF extrapolation using the HMI data [30]. Specifically, the total Lorentz force and torque are quantified by two surface integrals associated with three components of the photospheric magnetic field, and an optimization method is employed to minimize these two functions by modifying the magnetic field components within margins of measurement errors. Then Gaussian smoothing with FWHM of 8 arcsec is applied to all the three components of the magnetic field. Figure 1 compares the original HMI vector magnetogram with the preprocessed one as well as the final smoothed version. Using the vertical component B_z of this preprocessed and smoothed magnetogram, a potential field is extrapolated and used as the initial condition of the magnetic field in the MHD model.

In addition to the magnetic field, an initial plasma as the background atmosphere is also needed to start the MHD simulation. We used an isothermal plasma in hydrostatic equilibrium. It is stratified by solar gravity with a density ρ = 1 at the bottom and a uniform temperature of T = 1. Using this typical coronal plasma, we did not directly input the magnetic field into the model but multiplies it with a factor of 0.05 such that the maximum of magnetic field strength in normalized value is approximately 50~100 in the model. If using the original values of magnetic field, its strength (and the characteristic Alfvén speed) near the lower surface is too larger, which will be a too heavy burden on computation as the time step of our simulation is limited by the CFL condition. With the reduced magnetic field, we further modified the value of solar gravity to avoid an unrealistic large plasma β (and small Alfvén speed) in the corona. This is because if using the real number of the solar gravity (g_⊙ = 274 m s-2/gs=0.26), it results in a pressure scale height of H_p = 3.8, by which the plasma pressure and density decay with height much slower than the magnetic field. With the weak magnetic field strength we used, the plasma β will increase with height very fast to above 1, which is not realistic in the corona. To make the pressure (and density) decrease faster in the lower corona, we modified the gravity as

By this, we got a plasma with β < 1 mainly within z < 10 and the smallest value is 5 × 10⁻³.

Then we input the transverse field of the smoothed magnetogram to the model. This is done by modifying the transverse field on the bottom boundary incrementally using linear extrapolation from the potential field to the vector magetogram in time with a duration of 1 t_s until it matches the vector magnetogram. This will drive the coronal magnetic field to evolve away from the initial potential state, since the change of the transverse field will inject electric currents and thus Lorentz forces, which drive motions in the computational volume. Note that in this phase all other variables on the bottom boundary are simply fixed, thus the velocity remaining zero. This is somewhat un-physical since the Lorentz force will introduce non-zero flows on the bottom boundary, but it provides a simple and “safe” way (avoiding numerical instability) to bring the transverse magnetic field into the model. Once the magnetic field on the bottom surface is identical to that of the vector magnetogram, the system is then left to relax to equilibrium with all the variables (including the magnetic field) on the bottom boundary fixed. To avoid a too large velocity in this phase such that the system can relax faster, we set the kinetic viscosity coefficient as ν = 0.5Δx²/Δt, where Δx is the local grid spacing and Δt the local time step, determined by the CFL condition with the fastest magnetosonic speed. Actually this is the largest viscosity one can use with given grid size Δx and time step Δt, because the CFL condition for a purely diffusive equation with diffusion coefficient ν requires Δt ≤ 0.5Δx²/ν.

For the purpose of minimizing the numerical boundary influences introduced by the side and top boundaries of the computational volume, we used a sufficiently large box of (−32, −32, 0) < (x, y, z) < (32, 32, 64) embedding the field of view of the magnetogram of (−8.75, −8.25) < (x, y) < (8.75, 8.25). The full computational volume is resolve by a non-uniform block-structured grid with adaptive mesh refinement (AMR), in which the highest and lowest resolution are Δx = Δy = Δz = 1/16 (corresponding to 1 arcsec or 720 km, matching the resolution of the vector magnetogram) and 1/2, respectively. The AMR is controlled to resolve with the smallest grids the regions of strong magnetic gradients and current density. The magnetic field outside of the area of the magnetograms on the lower boundary is given as zero for the vertical component and simply fixed as the potential field for the transverse components. On the side and top boundaries, we fixed the plasma density, temperature, and velocity. The horizontal components of magnetic field are linearly extrapolated from the inner points, while the normal component is modified according to the divergence-free condition to avoid any numerical magnetic divergence accumulated on the boundaries.

In Figure 2, the curves colored in black show the evolution of the magnetic and kinetic energies integrated for the computational volume, and the residual of the magnetic field of two consecutive time steps which is defined as

where the indices k and k − 1 refer to the two consecutive time steps and i goes through all the mesh points. It can be seen that the magnetic energy increases sharply in a few time units², reaching ~1.4 of the potential field energy E₀ (here E0=1.2×1033 erg when scaled to the realistic value in the corona), and then keeps almost constant during the relaxing phase with bottom boundary fixed. Very similar, the kinetic energy first increases and later keeps on the level of 3 × 10⁻³ of E₀. The residual of the magnetic field increases in the first t_s as we continually modified the transverse field at the bottom boundary which drives quickly the evolution of the field in the corona. Then it decreases to below 10⁻⁵ with a time duration of 20t_s, which indicates that the magnetic field reaches a quasi-equilibrium state.

To make the field even closer to equilibrium, we carried out a “deeper” relaxation by running the model again, which is started with the relaxed magnetic field obtained at t = 20t_s and the initially hydrostatic plasma. Now we reduce the kinetic viscosity to ν = 0.05Δx²/Δt, i.e., an order of magnitude smaller than the previously used one, which will let the magnetic field relax further. Furthermore, the magnetic field at the boundary boundary is allowed to evolve in a self-consistent way by assuming the bottom boundary as a perfectly line-tying and fixed (i.e., v=0) surface of magnetic field lines. However, such a line-tying condition does not indicate that all magnetic field components on the boundary are fixed, because even though the velocity v is given as zero on the bottom boundary, it is not necessarily zero in the neighboring inner points. To self-consistently update the magnetic field, we solve the magnetic induction equation on the bottom boundary. Slightly different from the one in the main Equations (1), the induction equation at the bottom surface is given as

where we added a surface diffusion term defined by using a surface Laplace operator as ∇⊥2=∂2∂x2+∂2∂y2 with a small resistivity for numerical stability near the PIL ηstable=1×10-3e-Bz2, since the magnetic field often has the strongest gradient on the photosphere around the main PIL. The surface induction Equation (5) in the code is realized by second-order difference in space and first-order forward difference in time. Specifically, on the bottom boundary (we do not use ghost cell), we first compute v×B, and then use central difference in horizontal direction and one-sided difference (also 2nd order) in the vertical direction to compute the convection term ∇ × (v × B). The surface Laplace operator is also realized by central difference.

The curves colored in blue in Figure 2 show the evolution of different parameters during this relaxation phase³. Initially one can see a fast decrease of the magnetic energy because the magnetic field becomes more relaxed. As the viscosity is reduced significantly, the kinetic energy first increases, as driven by the residual Lorentz force of the magnetic field, to almost 5×10-3E0. Then as the magnetic field relaxed, the kinetic energy decreases fast to eventually less than 10-3E0, which is a very low level. The residual of magnetic field of two consecutive time steps also decreases to 10⁻⁵. These values show that the magnetic field reaches an excellent equilibrium.

Figure 3 shows the 3D magnetic field lines of the final relaxed MHD equilibrium. Note that the field lines are false-colored by the values of the force-free factor defined as α=J·B/B2, which indicates how much the field lines are non-potential. For a force-free field, this parameter is constant along any given field line. As can be seen, the magnetic field is close to a perfect force-free one since the color is nearly the same along any single field line. In the core of the configuration, the field lines are sheared significantly along the PIL, and thus have large values of α and current density. On the other hand, the overlying field is almost current-free or quasi-potential field with α~0, which plays the role of strapping field that confines the inner sheared core. Such a configuration is typical for eruption-productive ARs [e.g., 5, 31–34].

Based on the time sequence of vector magnetograms, it is straightforward to derive the surface velocity by employing the DARE4VM code developed by [19]. The differential affine velocity estimator (DAVE) was first developed for estimating velocities from line-of-sight magnetograms, and was then modified to directly incorporate horizontal magnetic fields to produce a differential affine velocity estimator for vector magnetograms (DAVE4VM). We use the SHARP data with cadence of 12 min and pixel size of 1 arcsec (by rebinning the original data with pixel size of 0.5 arcsec), and first of all, we fill the data gap using a linear interpolation in the time domain to generate a complete time series of 3 days from 00:00 UT on 2014 September 8 to 24:00 UT on September 10 with cadence of 12 min. Then we input the time series of vector magnetogram (rebinned as 1 arcsec per pixel) in the DAVE4VM code. A key parameter needed by the DAVE4VM code is the window size, and here we use 19 pixels, following [35] and [36]. After obtaining the surface velocity, we first reset those in the weak-field region (with total magnetic field strength below 100 G) as zero, because there are large errors and unresolved scales in these regions. Figure 4 shows evolution of the Poynting flux dE/dt, which is defined by

where S is the photospheric surface, and its time accumulation E_acc, as computed by using the surface flow (v_x, v_y, v_z) and the magnetic field. It can be seen that, except the data gap intervals⁴, the magnetic energy is continually injected in the corona through the photosphere, and in the 3 days, it gains ~5 × 10³² erg.

Before being input into the model, the flow data are also needed to be smoothed. We smoothed the time series of flow maps in both the time and space domains, with a Gaussian FWHM of 6 h (i.e., 30 time snapshots) and 8 arcsec, respectively, which is finally input to the data-driven model. Figure 5 shows 4 snapshots of the surface velocity after smoothing. The speed of the flow is generally a few hundreds of meters per second and the main feature is a clear and persistent rotation of the main sunspot. Note that during the 3 days the basic configuration of the photospheric magnetic flux distribution is rather similar with only somewhat dispersion. So the magnetic energy injection should come mainly from the transverse rotational flows. In addition to the rotational flow, we can see very evident diverging flow existing persistently near the boundary of the sunspot.

We input the surface velocity at the bottom boundary to drive the evolution of the model, by starting the simulation from the solution obtained from the time point of t = 58 in the relaxation phase (see Figure 2) as described in section 3. Here to save computing time, the cadence of the input flow maps, which is originally 12 min, was increased by 68.6 times when inputting into the MHD model. As a result, an unit of time in the simulation, t_s, corresponds to actually t_s × 68.6 = 7200 s, i.e., 2 h, in the HMI data. Compressing of the time in HMI data is justified by the fact that the speed of photospheric flows is often a few 100 m s⁻¹. So in our model settings, the evolution speed of the boundary field, even enhanced by a factor of 68.6, is still smaller by two orders of magnitude than the coronal Alfvén speed (on the order of 10³ km s⁻¹), and the quick reaction of the coronal field to the slow bottom changes should not be affected. The implementation of the bottom boundary conditions is the same as that for the deeper relaxation phase described in section 3. That is, on the bottom surface, we solved the Equation (5) to update all the three components of magnetic field with the flow field prescribed by those derived in section 4, while the plasma density and temperature are simply fixed. In the driven-evolution phase, the kinetic viscosity is also used as the smaller one ν = 0.05Δx²/Δt, which corresponds to a Reynolds number of 10 for the length of a grid cell Δx.

We show an approximately 2-day (a duration of 26t_s or 52 h in reality) evolution of the MHD system as an example, while the further evolution associated with an eruption and the physical mechanisms will be left for future study. Figure 6 presents the global energy evolutions. It can be seen the magnetic energy (the black line in the top panel) increases monotonously as driven by the surface flows. By the end of our simulation, it reaches approximately 1.8 times of the initial potential energy, and thus the total magnetic energy obtained in this driving phase is (1.8-1.4)×E0≈4.8×1032 erg. On the other hand, the kinetic energy (bottom panel of Figure 6) keeps below the level of 1×10-3E0 with a mild increase, which indicates that the system remains a stable, quasi-static evolution. In the top panel of Figure 6 we also plots the time integration of total Poynting flux (the blue line), using the magnetic field and the flow field on the bottom boundary of the simulation, which is the energy injected into the volume from the bottom boundary through the surface flow. If our boundary condition is accurately implemented, the energy injected from the bottom surface should match the magnetic energy obtained in the computational volume, since other energies are negligible if compared with the magnetic energy. As can be seen, the trend of magnetic energy evolution (the black line) matches rather well that of the energy input by the surface flow (the blue line), albeit a small numerical error that increases slowly the total magnetic energy, which is also seen in the relaxation phase as shown in Figure 2. It is worthy noting that the magnetic energy evolution is also in good agreement with the accumulated magnetic energy as derived from directly the observation data (the red line, which is the same value shown in Figure 4 but multiplied by a factor of 0.05² because the field input to the model is multiplied by 0.05) as calculated in section 4, which suggests that our data-driven model can reliably reproduce the magnetic energy injection into the corona through the photospheric motions. The small difference between the observation-derived and simulated energies might be due to the smoothing of the velocity since it filters out the small-scale flows that also contribute to the total Poynting flux.

Figure 7 (and its attached animation) shows the 3D magnetic field lines and their evolution in comparison with SDO/AIA image of coronal loops in the 94 Å wavelength which highlights the hot, core coronal loops in the AR. Overall, we can see a slow stressing of the field lines mainly through the sunspot rotational flow. It renders the core field lines to form a more and more coherent S shape, which resembles the observed sigmoid structure in the core of the AR. The increase of non-potentiality can also be seen from the increase of the force-free factor α in the core region. Figure 8 shows evolution of the surface magnetic flux distribution (see also the animation attached to Figure 7). The four snapshots are taken for the same times given in Figure 5, and thus comparing Figure 8 with Figure 5 shows the difference between the simulated magnetograms and the observed ones. In addition to the rotation of the main sunspot, an evident feature is the enhancement of the field strength along the PIL. This is owing to the divergence flow from the edge of the sunspot, which continually convect the magnetic flux to the PIL. Such pileup of magnetic flux near the PIL, however, is not seen in the observed magnetograms (Figure 5), which rather show field decaying. Such decaying of magnetic flux is likely due to the global turbulent diffusion of photospheric magnetic field by granular and supergranular convection [37] and other small-scale turbulence and flux cancellations in the photosphere, which is not being recovered by the DAVE4VM code and thus not reproduced in our simulation.

This paper is devoted to the description of a new approach of data-driven modeling of solar AR magnetic field evolution, in which, we have for the first time utilized directly the photospheric velocity field from DAVE4VM to drive a full-MHD model. To setup the initial conditions, we used a special MHD relaxation approach with sufficiently small kinetic viscosity to construct a true MHD equilibrium based on a single vector magnetogram. Then we derived the photospheric surface flows from a time series of observed magentograms based on the DAVE4VM method. The surface flows were finally inputted, again in time sequence, to the bottom boundary of the MHD model to self-consistently update the magnetic field at every time step, which is implemented by solving directly the magnetic induction equation at the bottom boundary using finite difference method.

We applied this data-driven model to study the magnetic field evolution of AR 12158 with SDO/HMI vector magnetograms. The initial MHD equilibrium is calculated using magnetogram observed for 00:00 UT on 2014 September 8, and a 2-day duration of the AR evolution is then simulated using the data-driven MHD model. Overall, the evolution is characterized by a slow stress of the field lines mainly through the rotational flow of the AR's leading sunspot, which makes the core field lines to forms an coherent S shape consistent with the sigmoid structure as seen in the AIA images. Such evolution proceeds in a quasi-static way since the kinetic energy of the system remains less than its magnetic energy by three orders of magnitude, while the magnetic energy increases monotonously as driven by the surface flow, and reaches approximately double of the initial potential energy by the end of the simulation. The magnetic energy obtained in the simulation during the surface driving phase matches closely the accumulated magnetic energy as calculated directly from the original vector magnetogram with the DAVE4VM derived flow field.

With the surface flow specified at the bottom boundary, the magnetic field can be updated self-consistently by solving the induction equation at the surface boundary. However, our simulation shows that discrepancy arises between the simulated magnetic field at the bottom surface with the original magnetograms. That is, in the simulation, magnetic flux are significantly enhanced along the PIL, owing to the divergence flow from the edge of the sunspot, whereas in the observed magnetograms such pileup of magnetic flux near the PIL is not seen, which rather shows field decaying. In the next step, we will try to use some ad-hoc flux cancellation near the PIL, for which a straightforward way is to increase the value of η_stable in Equation (5), such that the flux diffusion speed is comparable to that of the flux pileup. In the future, we will consider to include the global diffusion of photospheric magnetic field by granular and supergranular convection to simulate more realistically the magnetic field evolution. On the other hand, the discrepancy might result from errors in the vector magnetograms, since these errors can introduce spurious flows with the DAVE4VM code and thus influences our simulation. To elucidate this, we will test our model with error-free magnetograms from recent convective flux-emergence simulations [e.g., 38–40] as the ground-truth data.

Our ultimate purpose is to use the model to study how the coronal field, as driven by the slow photospheric motions, reaches a unstable state and runs into eruption. From this, we will be able to investigate in details the topology of the evolving magnetic field leading to the eruption, to see whether it forms a magnetic flux rope and becomes ideally unstable [e.g., 41, 42], or, a simply sheared arcade with a interally-formed current sheet to trigger flare reconnection [43], which would be helpful for resolving the long-standing debates on the triggering mechanism of solar eruptions [44].

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

CJ developed the numerical model, performed the result analysis, and wrote the draft. All authors participated in discussions and revisions on the manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest. The reviewer YG declared a past co-authorship with one of the authors CJ to the handling editor.