MHD Modeling of Solar Coronal Magnetic Evolution Driven by Photospheric Flow

It is well known that magnetic fields dominate the dynamics in the solar corona, and new generation of numerical modelling of the evolution of coronal magnetic fields, as featured with boundary conditions driven directly by observation data, are being developed. This paper describes a new approach of data-driven magnetohydrodynamic (MHD) simulation of solar active region (AR) magnetic field evolution, which is for the first time that a data-driven full-MHD model utilizes directly the photospheric velocity field from DAVE4VM. We constructed a well-established MHD equilibrium based on a single vector magnetogram by employing an MHD-relaxation approach with sufficiently small kinetic viscosity, and used this MHD equilibrium as the initial conditions for subsequent data-driven evolution. Then we derived the photospheric surface flows from a time series of observed magentograms based on the DAVE4VM method. The surface flows are finally inputted in time sequence to the bottom boundary of the MHD model to self-consistently update the magnetic field at every time step by solving directly the magnetic induction equation at the bottom boundary. We applied this data-driven model to study the magnetic field evolution of AR 12158 with SDO/HMI vector magnetograms. Our model reproduced a quasi-static stress of the field lines through mainly the rotational flow of the AR's leading sunspot, which makes the core field lines to form a coherent S shape consistent with the sigmoid structure as seen in the SDO/AIA images. The total magnetic energy obtained in the simulation matches closely the accumulated magnetic energy as calculated directly from the original vector magnetogram with the DAVE4VM derived flow field. Such a data-driven model will be used to study how the coronal field, as driven by the slow photospheric motions, reaches a unstable state and runs into eruptions.


INTRODUCTION
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 threedimensional (3D) coronal magnetic fields, numerical modelling 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 (Wiegelmann & Sakurai 2012), 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 modelling is becoming a viable tool to study coronal magnetic field evolution, which can selfconsistently describe both the quasi-static and the dynamic evolution phases (e.g., Wu et al. 2006;Feng et al. 2012;Cheung & DeRosa 2012;Leake et al. 2017;Inoue et al. 2018;Hayashi et al. 2019;Guo et al. 2019;Pomoell et al. 2019). Still, due to the limited constraint from observation, data-driven models are developed with very different settings from each other . For simplicity, some used the magneto-frictional model (Cheung & DeRosa 2012;Pomoell et al. 2019), 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 (Yang et al. 1986), although in some case an eruption can be reproduced but it evolves in a much slower rate than the realistic one (Pomoell et al. 2019), because the dynamic is strongly reduced by the frictional force. Some used the so-called zeroβ model (Inoue 2016;Inoue et al. 2018;Guo et al. 2019), 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 nonlinear 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 (Nakagawa 1980;Nakagawa et al. 1987;Wu & Wang 1987;Hayashi 2005;Wu et al. 2006) to specify the other variables according to the information of characteristics based on the wave-decomposition principle of the full MHD system (e.g., . 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 (Schuck 2008) to drive the MHD model. The differential affine velocity estimator (DAVE) was first developed for estimating velocities from line-ofsight 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 (Priest & Forbes 2002). 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 projectedcharacteristic 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 −1 , which is well documented in the literature (Vemareddy et al. 2016;Zhao et al. 2016;Zhou et al. 2016). The AR is well isolated from neighboring ARs, thus is suited for our modelling 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 nonpotential 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 five days (Vemareddy et al. 2016). 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., Vemareddy et al. 2016;Zhao et al. 2016;Zhou et al. 2016), while some nonlinear force-free field extrapolations appear not to support this (Duan et al. 2017;Lee & Magara 2018). 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.

MHD EQUATIONS
We numerically solve the full MHD equations in a 3D Cartesian geometry by an advanced conservation element and solution element (CESE) method (Jiang et al. 2010). Before describing the model equations in the code, it is necessary to specify the quantities used for nondimensionalization. Here we use typical values at the base of the corona for non-dimensionalization as given in Table 1 Table 1. Parameters used for non-dimensionalization. n is a typical value of electron number density in the corona given by n = 1 × 10 9 cm −3 and m is the mean atomic mass.
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 −ν ρ (ρ − ρ 0 ) to the continuity equation, where ρ 0 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 everdecreasing 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 v A = 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 sim-plicity, such that the energy equation describes an isothermal process. 1 The kinetic viscosity ν will be given with different values when needed, which is described in the following sections.

CONSTRUCTION OF AN INITIAL MHD EQUILIBRIUM
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 smallscale 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., Wiegelmann 2004), 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 ). The preprocessing is done using a code developed by Jiang & Feng (2014), which is originally designed for NLFFF extrapolation using the HMI data (Jiang & Feng 2013). 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 /g s = 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 −3 . 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 nonzero 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 2 /∆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 2 /ν.
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 divergencefree 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 2 , reaching ∼ 1.4 of the potential field energy E 0 (here E 0 = 1.2 × 10 33 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 −3 of E 0 . 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 −5 with a time duration of 20 t 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 = 20 t s and the initially hydrostatic plasma. Now we reduce the kinetic viscosity to ν = 0.05∆x 2 /∆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 selfconsistent 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 ∂x 2 + ∂ 2 ∂y 2 with a small resistivity for numerical stability near the PIL η stable = 1 × 10 −3 e −B 2 z , 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 secondorder 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 The curves colored in blue in Figure 2 show the evolution of different parameters during this relaxation phase 3 . 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 −3 E 0 . Then as the magnetic field relaxed, the kinetic energy decreases fast to eventually less than 10 −3 E 0 , which is a very low level. The residual of magnetic field of two consecutive time steps also decreases to 10 −5 . 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 falsecolored by the values of the force-free factor defined as α = J · B/B 2 , 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 currentfree 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., Schrijver et al. 2008;Sun et al. 2012;Toriumi & Wang 2019;Duan et al. 2019).

DERIVE THE SURFACE FLOW FIELD
Based on the time sequence of vector magnetograms, it is straightforward to derive the surface velocity by employing the DARE4VM code developed by Schuck (2008). 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 minutes 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 three days from 00:00 UT on 2014 September 8 to 24:00 UT on September 10 with cadence of 12 minutes. 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 Liu et al. (2012) and Liu et al. (2014). 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 4 , the magnetic energy is continually injected in the corona through the photosphere, and in the three days, it gains ∼ 5 × 10 32 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 hours (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 three 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.

DATA-DRIVEN SIMULATION
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 hours, 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 −1 . 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 3 km s −1 ), 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 2 /∆t, which corresponds to a Reynolds number of 10 for the length of a grid cell ∆x.
We show an approximately two-day (a duration of 26 t s or 52 hours 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) × E 0 ≈ 4.8 × 10 32 erg. On the other hand, the kinetic energy (bottom panel of Figure 6) keeps below the level of 1 × 10 −3 E 0 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 2 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 (Wang et al. 1989) 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 modelling 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.

CONCLUSIONS
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 two-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 poten-tial 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., Chen et al. 2017;Cheung et al. 2019;Toriumi & Hotta 2019) 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., Kliem & Török 2006;Aulanier et al. 2010), or, a simply sheared arcade with a interally-formed current sheet to trigger flare reconnection (Moore et al. 2001), which would be helpful for resolving the long-standing debates on the triggering mechanism of solar eruptions (Chen 2011).