{'fraction_non_alphanumeric': 0.026728032046364952, 'fraction_numerical': 0.01585272308872411, 'mean_word_length': 4.941462426574843, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Data-driven modeling of solar coronal magnetic field evolution and eruptions

Magnetic fields play a fundamental role in the structure and dynamics of the solar corona. As they are driven by their footpoint motions on the solar surface, which transport energy from the interior of the Sun into its atmosphere, the coronal magnetic fields are stressed continuously with buildup of magnetic nonpotentiality in the form of topology complexity (magnetic helicity) and local electric currents (magnetic free energy). The accumulated nonpotentiality is often released explosively by solar eruptions, manifested as solar flares and coronal mass ejections, during which magnetic energy is converted into mainly kinetic, thermal, and nonthermal energy of the plasma, which can cause adverse space weather. To reveal the physical mechanisms underlying solar eruptions, it is vital to know the three-dimensional (3D) structure and evolution of the coronal magnetic fields. Because of a lack of direct measurements, the 3D coronal magnetic fields are commonly studied using numerical modeling, whereas traditional models mostly aim for a static extrapolation of the coronal field from the observable photospheric magnetic field data. Over the last decade, dynamic models that are driven directly by observation magnetograms have been developed and applied successfully to study solar coronal magnetic field evolution as well as its eruption, which offers a novel avenue for understanding their underlying magnetic topology and mechanism. In this paper, we review the basic methodology of the data-driven coronal models, state-of-the-art developments, their typical applications, and new physics that have been derived using these models. Finally, we provide an outlook for future developments and applications of the data-driven models.


INTRODUCTION
The solar corona, the upper atmosphere of the Sun, consists of highly magnetized, hot plasma with temperatures of millions of degrees. 1 Observations using soft X-ray and extreme ultraviolet (EUV) emissions of this plasma indicate that the magnetic field in the corona has a complex and evolving structure that is continuously driven by the ceaseless motions at the solar surface, the photosphere, by which energy is transported from the interior of the Sun to the corona. During this process, the coronal magnetic field is stressed, electric current is built up, and magnetic free energy is accumulated continuously. From time to time, magnetic explosions, known as solar eruptions, are produced as a result of the rapid release of this magnetic energy into thermal and kinetic energy. 2,3 The most violent eruptions manifest as major solar flares and coronal mass ejections (CMEs). They unleash a huge amount of magnetic flux (on the order of 10 14 Wb or 10 22 Mx, the same order of magnetic flux contained in a typical solar active region) and plasma (with a mass of 10 11 $ 10 13 kg, similar to the mass of a typical solar filament or prominence) into the heliosphere 4 and drive adverse space weather events that could severely affect human activities in modern society.
It has long been recognized that magnetic fields play a fundamental role in almost all dynamics in the solar corona, such as formation of coronal loops and prominences (or filaments), initiation of solar eruptions, as well as modulation of the global structure of solar wind. 5 However, it remains very difficult to make a direct measurement of the coronal magnetic fields from emissions of the extremely tenuous plasma in the corona by traditional spectro-polarimetric methods (using the Zeeman and Hanle effects). There are attempts made to measure coronal fields using infrared wave bands 6 and radio emissions [7][8][9] as well as coronal seismology 10 (or magnetoseismology 11,12 ), but these can only give fragmentary and occasional data. Another way has been proposed to measure coronal fields, using magnetically induced transitions of EUV lines (e.g., Fe X and Fe XI), but the accuracy of determination of field strength would be of the order of the magnitude. 13,14 Up to the present, routine measurement of the Sun's magnetic field that we can rely on is restricted to only a single layer of the solar atmosphere, i.e., the photosphere (the photosphere still has a finite thickness, but it is extremely small compared with the length scales in the corona).
The concept of numerically modeling the three-dimensional (3D) coronal magnetic field with the bottom boundary constrained by the observed photospheric magnetograms dates back to over half a century ago. 15 The earliest model used a lowest-order assumption that the coronal field is current free (namely, a potential field), whereas almost all modern codes for coronal magnetic field reconstruction (or extrapolation) are based on the more realistic assumption that the Lorentz force vanishes in the corona (i.e., the magnetic pressure force and magnetic tension force balance each other). This is because in most parts of the low corona (with a height of less than $400 Mm from the coronal base), particularly in the strongly magnetized active regions (ARs), the plasma b (the ratio of plasma thermal pressure to magnetic pressure) is low (with typical values of 10 À4 to 10 À2 ), and all other forces can be neglected compared with the magnetic pressure gradient force (or magnetic tension force). Over the past several decades, a variety of nonlinear force-free field (NLFFF) extrapolation codes that use vector magnetograms as the key input have been developed and have greatly improved our knowledge of the 3D structure and topology of the coronal magnetic fields, which is extensively discussed in a few review papers. [16][17][18][19] The NLFFF models have major limitations, and one of them is due to the basic assumption that the coronal magnetic field always remains in a static equilibrium. Such a Lorentz force balance is strongly violated during solar eruptions, in which the magnetic field evolves rapidly with Alfvénic speed as driven by the unbalanced Lorentz force. The interaction of the magnetic field with plasma is totally omitted in the NLFFF and other similar static extrapolation models, and thus they cannot be used to study dynamics associated with the temporal evolution of the corona. Another limitation of some NLFFF models, especially of those employing an iteration approach and assuming a potential field as the initial guess, would be that it may be difficult to reach the solution when the actual coronal fields are highly deviated from the potential field (such as those including magnetic flux ropes).
The macroscopic behavior of coronal dynamic evolution is described selfconsistently by magnetohydrodynamics (MHD) equations, which characterize the nonlinear interaction of the magnetic field with plasma. 20,21 In recent years, time-dependent models based on MHD (or its simplifications) and using continuously observed magnetograms as input have arisen to reproduce the dynamics of the corona. Such data-driven models are a significant step forward in coronal magnetic field modeling compared with traditional methods, providing a novel avenue to study coronal magnetic evolution, including slow, quasi-static coronal evolution as driven by photospheric motions, the rapidly evolving field in solar eruptions, and the self-consistent transition from a quasi-static stable equilibrium to fast eruption. In this paper, we review the basic methodology of data-driven coronal models, state-of-the-art developments, their typical applications, and new physics that have been derived using these models. This review is focused on low coronal evolution and mainly on the scale of ARs with emphasis on solar eruptions.

The magnetofrictional model
The magnetofrictional model was originally proposed for NLFFF reconstruction. 51 It is the computationally cheapest way (and, thus, the most frequently used one) to carry out data-driven simulations by assuming that the coronal magnetic field evolves slowly as a continuous series of force-free equilibria. To derive this model, an artificial frictional force v f rv (where v f is the frictional coefficient, r the plasma density, and v the velocity) is included in the MHD equation of motion, where p, g, j, and B are the thermal pressure, the solar gravity, the electric current density, and the magnetic field, respectively. Then, by neglecting the thermal pressure, the gravitational force, and the inertial effect, this equation can be reduced to expression of the so-called magnetofrictional velocity: The magnetofrictional velocity acts to ensure that the magnetic field remains close to a force-free equilibrium as the field is perturbed via boundary motions. A pseudo-density is assumed to be r = B 2 (thus, a uniform pseudo Alfvén speed; the pseudo density is usually absorbed into the frictional coefficient), and the only equation that needs to be solved is the induction equation (Faraday's law), where E is the electric field, and hj is the resistivity term which can be omitted for an ideal process. To fulfill the solenoidal requirement of B, the vector potential A (with a Coulomb gauge; i.e., V,A = 0) can be used as the primary variable (thus, B = V 3 A), and the induction equation is then written as 23,26,32 vA vt = À E = v 3 B À hj: A central difference scheme with a staggered grid 52 can be used to calculate the spatial derivatives, and then a forward difference in time is used to advance the solution. The time step should be determined by the largest Alfvén speed, complying with the Courant-Friedrichs-Lewy (CFL) condition. There are different choices for the frictional coefficient, which determines the relaxation speed of the system (some use a constant value 26,32 and others use a height-dependent value 22,23,27,28 ) to optimize the match between simulations and observations. The initial value of the magnetic field can be obtained from static extrapolations (e.g., a potential field or more sophisticated NLFFF if needed).
Compared with other models, the magnetofrictional model is a convenient way to follow the quasi-static coronal evolution, but it should be noted that this model, in principle, cannot be used to simulate eruptions whose dynamics are far away from a force-free state. Another serious issue of the magnetofrictional model, as recently pointed out by Yeates, 53 is that it might lead to breakdown of magnetic flux conservation even when using an ideal induction equation.

Zero-b MHD model
The zero-b MHD model is also a simplified version of the MHD model. It omits the non-Lorentz forces of gas pressure gradient and solar gravity in the equation of motion, which thus reads as where a viscosity term is often included for the purpose of fast relaxation (and numerical stabilization). In some works, a pseudo-density similar to that in the magnetofrictional model (for example, r = B 2 or r = B ) is used; 37 with an initial ad hoc setting of the density like r 0 fB 3=2 0 (where B 0 is the initial magnetic field strength), aimed at mimicking the profile of the coronal Alfvén speed, 56 or a more realistic one that employs the stratified density profile from the photosphere to the corona; 35 thus, the observed photospheric magnetic field can be handled at the bottom boundary directly. The zero-b model can be solved using existing full MHD codes; for example, the data-driven model developed by Guo et al. 35 is built on an open-sourced parallelized, adaptive mesh refinement (AMR) code for MHD, MPI-AMRVAC. 57,58 By discarding the thermal pressure (and, thus, the energy equation), the zero-b model can be solved more efficiently than the full MHD model, and it avoids the particular numerical problem of negative pressure potentially arising in the very low-b region. But the caveat is that the zero-b assumption is reasonable only in domains where the plasma b is extremely low, whereas the real b in the corona is quite variable; 59 for example, there are also places with considerably high b in the corona because of the prevalence of magnetic null points. Relevant to this limitation is that, for studying eruptions, the zero-b assumption might fail when there is a fast reconnection in the field, in which the thermal pressure could play an important role in dynamics in the weak-field region of magnetic field dissipation.

The full MHD model
The full MHD model solves all equations, including the continuum equation, the momentum equation, the induction equation, as well as the energy equation, and thus has none of the limitations mentioned for the magnetofrictional and zero-b models, but, of course, the computational cost will be higher. The full MHD equations as used by different authors differ mainly in the last equation (i.e., the energy equation), which can take various forms, using different assumptions for particular purpose. For instance, an adiabatic process is often considered when focusing on the evolution of the magnetic field and its interaction with plasma; e.g., using a temperature equation, 44,60 vT vt where T is the temperature and g the polytropic index. By letting g = 1, it reduces to the isothermal process. 61 More comprehensive thermodynamics could be included, such as the radiation, thermal conduction, and (parameterized) coronal heating processes, 42,47,62,63 to build data-driven radiative MHD model. 64,65 For localized simulation of AR scales, the initial conditions often consist of an isothermal 44,60 or temperature-stratified 42,47,66 atmosphere in a hydrostatic state and a potential or near force-free magnetic field. Aimed at mimicking the highly non-uniform distribution of plasma in the corona, some models 48 also use a pre-initial density proportional to local magnetic field strength and let the MHD system evolve to a new equilibrium, which is then used as the initial state of subsequent data-driven simulation. For a global simulation that includes solar wind, the initial conditions are usually obtained by relaxing the Parker's spherically symmetric solution of solar wind along with a global potential field to a steady state. Although the classical finite difference method can be used, 42,47 numerical integration of the MHD equations is more commonly implemented using the modern, Godnuov-type finite volume scheme; for example, the total variation diminishing (TVD) Lax-Friedrichs formulation 67 and the monotonic upstream scheme for conservation laws 68 based on characteristic decomposition and Riemann solver. A recently developed scheme, the space-time conservation element and solution element (CESE) method, [69][70][71] has also been used efficiently in data-driven MHD models with features of massive parallel computing and AMR. 41,44,45 THE DATA-DRIVEN BOUNDARY CONDITIONS The unique feature that distinguishes a data-driven coronal model from others is incorporation of the observables of a time sequence directly into the model. This is realized by specifying time-dependent boundary conditions at the bottom surface using the continuously observed data as provided there. Ideally, all independent variables that are required to be solved in the controlling equations of the specific models should be provided self-consistently for the bottom boundary by observations. 72 However, because of the limited observations, these boundary conditions are underspecified, and, thus, different implementations use different assumptions, which currently can be classified into mainly three types. The first choice, B driven, is to give the observed photospheric magnetic field data directly as the boundary values. The second one, V driven, is to determine the plasma flow at the photosphere and to drive the simulated magnetic field through the Review ideal induction equation, v t B = V 3 ðv 3 BÞ, or other ways. The third choice, E driven, is to determine the electric field E and to drive the magnetic field as v t B = À V 3 E. A combination of different choices also exists; for example, by specifying simultaneously the observed magnetic field and the derived velocity at the bottom boundary. A detailed description of the implementation of these boundary conditions is given in the supplemental information.

DATA-DRIVEN MODELING OF THE SLOW EVOLUTION OF CORONAL MAGNETIC FIELDS
In a period without eruptions, the coronal magnetic fields evolve in a quasistatic way, driven by the slow motions at the photosphere. Although such an evolution could be approximately followed by a time series of independent NLFFFs 73-75 or MHD equilibria, 40,49,76,77 the data-driven model produces a continual evolution that allows a "memory" of magnetic connectivity and buildup of electric currents and free magnetic energy. Wu et al. 48,78 studied the realistic evolution of ARs using a data-driven full MHD model for the first time. Their model is driven by the line-of-sight (LOS) magnetograms from the Michelson Doppler Imager (MDI) onboard Solar and Heliospheric Observatory (SOHO) and includes the effects of differential rotation and meridional flow, and concluded that, for the studied AR, the surface transverse flow plays a dominant role in magnetic energy flux through the photosphere.
Because the full MHD simulation is time consuming, the computationally faster magnetofrictional models, which utilize LOS or vector magnetograms, are used widely to follow the long-term quasi-static evolution of ARs; e.g., a few days during which ARs build up their nonpotentiality. Some of the works are focused on the magnetic energy injection from the photosphere into the upper atmosphere, which can shed light on the coronal heating problem. 79 Mackay et al. 26 modeled the decaying phase of AR NOAA 8005 using SOHO/ MDI magnetograms, quantified the free energy accumulation as driven by the magnetic flux dispersal, and concluded that the injected free energy is sufficient to explain the radiative losses at coronal temperatures within the AR. Using data from the Helioseismic and Magnetic Imager (HMI) onboard Solar Dynamics Observatory (SDO)_and data from Sunrise balloon-borne observatory, Chitta et al. 80 applied a similar model to study the spatial and temporal variations of the energy dissipation rate and energy flux in a quiet region. They found that the energy deposited in the solar atmosphere is concentrated mostly near the photosphere (within a height of only 2 Mm) and that there is no sufficient energy flow upward to balance radiative and conductive losses in the corona. The datadriven magnetofrictional model has also been extended to the global corona to study the non-potential energy evolution of the global magnetic field 30 and to reproduce the largescale magnetic topology with full-disk EUV images ( Figure 1) for the middle corona 81 over a few Carrington rotations.
Other studies using the magnetofrictional models are devoted to formation and evolution of the eruptive structure, in particular, magnetic flux rope (MFR), which is commonly believed to be an important candidate for the pre-eruptive magnetic structures. [83][84][85][86][87][88][89] Cheung and DeRosa 23 performed simulations of the flux-emerging AR NOAA 11158 and found that, when continuous twisting is imposed at the bottom boundary (through an assumed form of electric field), the AR can form a series of flux rope ejections from the central polarity inversion line (PIL). Using a proxy emissivity based on the LOS integration of the field-line averaged square of current density, they were able to produce synthesis images with a visual texture similar to EUV images of coronal loops ( Figure 2). Cheung et al. 22 modeled a recurrent MFR formation during the emergence of a currentcarrying magnetic field that produced homologous helical jets. Gibb et al. 24 modeled how a sigmoidal MFR is built up in the decaying phase of AR NOAA 10977 through continual flux cancellation. 90 Yardley et al. 32 simulated nearly the full life of AR NOAA 11437 (5 days) from its initial emergence to decay phases, in which MFRs were also found to form and erupt. By parameterized tests of different non-inductive electric field entered at the boundary for simulation of AR NOAA 11504, Pomoell et al. 27 showed that MFRs can form and erupt only with an enhanced injection of magnetic helicity. They also found that the coronal evolution appears to be significantly different with different helicity injections despite relatively similar total energy injections because the helicity quantifies the overall degree of complexity of the field structure. Price et al. 28 used the same data-driven approach to study formation of an unstable MFR in the rapid and complex flux-emerging AR NOAA 12673, which produced the largest flare in solar cycle 24. An interesting finding is that the evolution of magnetic helicity in the simulation suggests a possible threshold for eruptions in the ratio of current-carrying helicity to relative helicity. 91,92 Price et al. 29 further simulated an MFR formed in AR NOAA 12473 and found it became torus unstable 93 before the eruption. Kilpua et al. 25 analyzed the details of an MFR eruption in AR NOAA 12158 and attempted to compare its structure (but still in the source region) with in situ observations at 1 AU.
Although many studies using the data-driven magnetofrictional model also claimed simulation (or even reproduction) of eruptions, the magnetofrictional model cannot self-consistently model the eruptive dynamics because of the inherent limitations of the model; i.e., the assumption that the coronal field evolves quasi-statically. First, the strong frictional force in the model can quickly slow down any unstable structure that tends to erupt, especially considering that there is no inertial effect. Thus, the erupting structure moves very slowly. For example, as shown in the simulation of a torus-unstable MFR by Price et al., 29 the MFR takes over 4 days to ascend a height of about 100 Mm, which otherwise should only take a few minutes in a realistic eruption. Second, the fast reconnection, 95 which is responsible for the impulsive release of free magnetic energy into kinetic energy, which powers the eruption, cannot happen in the magnetofrictional model. So these simulations ignored the important role of fast reconnection in driving the eruption that is proposed in many theoretical models of eruption. [96][97][98][99] Therefore, it is more appropriate to use the pre-eruptive field as simulated in the magnetofrictional model as an initial condition to a more self-consistent full MHD model for the purpose of simulating eruptions. 100

SIMULATION OF SOLAR ERUPTIONS
The MHD simulation initialized with unstable (or close to unstable) magnetic field data that are extrapolated for a time shortly prior to an eruption proves to be a very useful approach to realistically reproduce the subsequent eruptive process, which has been shown by Jiang et al. 101 when studying the eruption of a sigmoidal MFR in AR NOAA 11283. Strictly speaking, such a simulation cannot be considered a data-driven model because only a single set of magnetograms is used (thus, data constrained would be more appropriate), but we regard it as a quasi-data-driven model by considering that the initial condition is constructed using data, and it is the initial instability that drives the evolution. Because, with just a single snapshot of a magnetogram, such a type of model is only valid for a short period (e.g., a few minutes, focusing only on the eruption process), during which the photospheric driving effect (e.g., the surface motion and flux emergence and cancellation) can be neglected.
Analysis of the pre-flare NLFFF extrapolated from the SDO/HMI vector magnetogram based on their CESE-MHD-NLFFF code, 54,102,103 Jiang et al. 101 show that AR NOAA 11283 has a spine-fan topology of magnetic null point, linking multiple polarities and an MFR, which is close to torus instability (TI), 93 is embedded in the complex topology. When entered in the MHD model (although with a fixed bottom boundary), the field erupts, driven by the TI, and is accelerated by the breakouttype reconnection 96 at the null. The initiation and subsequent evolution of the eruptive flux rope resemble remarkably in morphology the corresponding EUV features, such as the sigmoid, the filament ejection, as well as quasi-circular flare ribbons, as observed by the Atmospheric Imaging Assembly (AIA 104 ) on board SDO ( Figure 3). This result has two-fold significance. On one hand, it evidences the agreement of the extrapolated field with the realistic corona prior to eruption; on the other hand, it shows the reliability of the MHD model in characterizing the dynamics during eruptions. Because previous simulations of solar eruptions before Jiang et al. 101 were commonly carried out with idealized or hypothetical magnetic field configurations, the combination of NLFFF and MHD models provides a viable way to reproduce realistic eruption process with complexity that is beyond the scope of idealized or theoretical models.
A similar way of simulating realistic eruptive events has been employed by many modelers with different pre-flare coronal field reconstructions (on which the performance of such kind of simulation largely depends), different MHD codes, as well as different boundary conditions. [35][36][37]39,73,100,[105][106][107][108][109][110][111][112] For example, Kliem et al. 110 studied the first eruptive event witnessed by SDO using an unstable NLFFF constructed with the flux rope insertion technique 113,114 to initialize a zerob MHD model, which yields good agreement with the strongly inclined (non-radial) rise path of the filament. The same approach is used to reveal the process of magnetic twist release during a filament eruption. 112 Similarly, Guo et al. 34 studied another filament eruption using their zero-b MHD model with the initial condition of the magnetic field determined by an advanced flux rope insertion method that is based on regularized Biot-Savart laws. The flux rope insertion method can reconstruct the MFR reasonably well, even in a weak-field region for which the traditional NLFFF extrapolation codes often fail because of the large uncertainty of measurements in the vector magnetograms, but fine-tuning of the free parameters is often required to reach a solution consistent with the observed filament. [115][116][117][118] There are also a few studies 111,119-121 that use a non-force-free coronal field extrapolation 122-124 as the initial conditions of the MHD simulation, in which the eruption is directly triggered and likely driven by the initial unbalanced Lorentz force.
Some works also reconstruct a coronal magnetic field that is slightly farther away from being unstable and drive it to eruption using ad hoc ways that are not specified by the observed data. Inoue et al. 36,106,107 investigated the eruption mechanism of eruptive flares in the well-known AR NOAA 11158. 75 They first extrapolated NLFFF using SDO/HMI vector magnetogram observed for around 2 h before the flares, and found the NLFFF is stable in their zero-b MHD simulation. Then they applied enhanced anomalous resistivity to trigger tether-cutting recon-nection 125 in the sheared arcades of the AR core, which eventually builds up a torus-unstable MFR to eruption. It is stressed 36 that, after the initiation phase, the eruption is mainly driven by the ongoing reconnection rather than TI, which explains why eruptions can still succeed even in a background field with a saddle-like profile of the decay index (i.e., the controlling parameter of TI with a value of around 1.5 as its threshold 84,[126][127][128][129], first with a TI-unstable and then a TI-stable area (Figure 4). With vector magnetograms provided by the Solar Optical Telescope (SOT 130 ) onboard Hinode, 131 Amari et al. 73   MFR is found to be formed gradually. Then the NLFFF based on the last available vector magnetogram, which is taken around 6 h before the flare is entered into the MHD model. Like the work of Inoue et al., 106 the input NLFFF model is still stable and shows no eruption in the MHD model. Further stress of the coronal field through artificially specified converging flows on the boundary conditions, which mimics flux cancellation, is used to drive the MFR to approach an unstable state till its eruption. Using the same approach, Amari et al. 135 simulated the failed eruption of an MFR in the super AR NOAA 12192 (the largest AR in the last two solar cycles 136 ), in which the so-called "magnetic cage" (i.e., the overlying strapping field lines) was found to play a key role, resulting in the unusual confined flares of X-classes in that AR. 137 The eruption of AR NOAA 10930 was also studied by Muhamad et al., 138 who first extrapolated a stable NLFFF model for the pre-flare field and then let a small bipole emerge into the large-scale field with different orientations relative to the main PIL. They confirmed that the bipole field with a so-called opposite polarity and reversed shear structures can effectively trigger flare eruption, which was originally suggested by Kusano et al. 139 With realistic simulations that are constrained or driven by observed data, significant insight into complex eruption events can be gained. By an advanced simulation based on their previous work, 101 Jiang et al. 140 demonstrated a three-stage magnetic reconnection scenario of the sigmoid eruption in AR NOAA 11283. These three reconnection episodes occurred successively in different locations in the corona and led to a complex process of CME formation. The initial sigmoidal flux rope breaks one of its legs, and quickly gives birth to a new tornado-like magnetic structure that is highly twisted and has multiple connections to the Sun at remote places from the original flare site. The final weak reconnection triggered in the newly formed current sheet associated with the relatively weak field that connects the remote polarity might provide a viable explanation for the EUV later phase of the flare. 141,142 Inoue et al. 37,108 simulated the X9.3 flare in AR NOAA 12673 and revealed a complex eruption process in which a series of small flux ropes straddled along the same PIL prior to the flare reconnect to form a large, highly twisted MFR that might trigger its eruption by kink instability. 143 The same event was also simulated by Jiang et al., 109 who focused on accurately char-acterizing the evolving magnetic topology during the eruption and were able to reproduce the spatial location and the temporal separation of the complex flare ribbons and the dynamic boundary of the erupting MFR's feet. The temporal profile of the total reconnection flux is comparable with the soft X-ray light curve. A recent ultra-high resolution, fully 3D MHD simulation 144 (reaching a spatial scale of 45 km) even reproduced the successive formation of mini flux ropes (i.e., plasmoids in 2D) in the reconnection of a confined flare that matches the high resolution from New Vacuum Solar Telescope (NVST) 145 and SDO, as shown in Figure 5.
Although with such success, the combination of static reconstruction for preeruptive field and time-dependent simulation of its subsequent eruption cannot self-consistently show how the pre-eruptive field is formed and how the eruption is triggered (but such simulations can be used to test the different scenarios). One of the issues is that many NLFFF models contain a numerical residual of Lorentz forces, which is often non-negligible, as indicated by the misalignment of the current density and magnetic field; thus, the initial magnetic field already deviates from a force-balanced state because of the numerical residual forces. For example, the current-weighted sine of the angle between J and B for many NLFFF models is typically in the range of 0.2-0.4. 147-150 Such residual force can instantly induce plasma motions in a low-b and highly tenuous plasma environment. This motion provides a perturbation to trigger instability and, thus, eruption of the magnetic field if it is unstable or not far from unstable. But there is also the possibility that the large residual force may overtake the instability, and, thus, the true initiation mechanism of eruption may not be identified. So, it is more promising when the NLFFF can be relaxed to a well-established equilibrium by cleaning the residual forces and to use this equilibrium as the initial conditions of MHD simulation. Although there are some studies 73,106 beginning simulation from a stable NLFFF, they need to artificially drive the field to evolve toward eruption by ad hoc boundary flows or enhanced resistivity. In this case, it is generally not guaranteed that the resulted evolving magnetic field at the bottom boundary will match the evolution of the observed magnetogram. In other words, the selfconsistent evolution of the coronal field from a stable equilibrium to eruption matching the variation in the photospheric field is still out of reach.  36 ). The dashed line denotes the critical decay index of n = 1.5 (i.e., the threshold of TI), and the light-blue and blue masks mark regions n > 1.5 and n < 1.5, respectively. The flux rope can erupt successfully even with the decay index of a saddle-like profile. Image reproduced from Inoue et al. 36

DATA-DRIVEN SIMULATION OF QUASI-STATIC EVOLUTION TO ERUPTION
A self-consistent, data-driven model of solar eruptions should be able to simulate the coronal magnetic field following a long-duration, quasi-static evolution all the way to its fast eruption. Jiang et al. 44 developed such an MHD model (the socalled driven ARE (DARE)-MHD model 44,151,152 ) and successfully simulated a flux emergence event of over 3 days that finally leads to an eruption (of an M-class flare with a CME) in the topology-complex AR NOAA 11283. The simulation is started with a potential field extrapolation from the vertical component of the magnetogram taken for the initial time; thus, one can see how the pre-flare magnetic free energy is accumulated continually from the very beginning, driven by the flux injection with a strong shearing motion. The simulated field morphology well resembles the AIA observations for the whole process. The timing of the simulated eruption onset matches the observed flare time within an offset of less than 2 h (by comparing it with the whole simulation of 3 days), which supports the model capturing the key transition of dynamics from pre-eruption to eruption. By a thorough analysis of the magnetic topology evolution, Jiang et al. 44 concluded that the flux emergence leads to formation of a jet-like configuration that is favorable for a self-amplifying reconnection between the newly emerged non-potential arcade and the pre-existing open flux, which is subsequently triggered and results in the eruption (Figure 6). This scenario is actually similar to the breakout model. 96 In this event, no MFR formed, but while studying another event, the DARE-MHD simulation supports the high-resolution observations by the Goode Solar Telescope 153 and SDO of flux emergence of a small-scale MFR and its subsequent eruption. 154 Similar simulations of the large-scale MFR formation and eruption have been carried out for CME initiation in AR NOAA 12371. 46,151 The DARE-MHD model has also been used to study confined flares. For instance, Jiang et al. 152 simulated a 2-day magnetic evolution that lead to the X3.1 confined flare in the super AR NOAA 12192. They showed that a thin layer of intense current is formed gradually in the AR's core. The current layer was successively enhanced, driven by the photospheric shearing motion, until it became so thin (as a current sheet) that a tether-cutting reconnection 125 between the sheared magnetic arcades sets in, which is responsible for the flaring process, and after the reconnection, the current is substantially dissipated. Comparison with the AIA observations shows that the footpoints of the newly reconnected field lines reproduced almost exactly the location of the flare ribbons, and the morphology of the reconnecting field lines and the simulated EUV image (using a similar method as proposed by Cheung and DeRosa 23 ) well resembles the flaring loops. According to the simulation, this flare failed to erupt as a CME because the consequence of the reconnection is still a sheared arcade without forming an escaping MFR because the pre-flare field is only weakly twisted, 155 but of course the strong overlying field also plays an important role. 135 Recently, efforts have been made to compare the simulated eruption with the observed features more quantitively. For example, Guo et al. 35 simulated a shorttime evolution of about 1 h before a small flare (C4.7) to its onset in AR NOAA 11123 using their data-driven zero-b MHD model and reproduced an MFR erup-tion with a morphology highly consistent with the corresponding AIA observations of the flare-associated filament eruption (Figure 7). A detailed analysis of the kinematics of the rise of the MFR showed that the simulated eruption distinctly demonstrates a two-stage evolution similar to the slow-rise and fast-rise phases of many observed eruptions; 156 by computing the decay index, they suggested that the slow rise to fast eruption is likely due to TI. Using the same model, Zhong et al. 39 showed the capability of simulating a failed eruption that occurred on January 30, 2015. Similarly, their zero-b model is initialized with an NLFFF extrapolated for a short time of half an hour before flare onset and successfully recreated the initial rise of an MFR (likely because of kink instability) and then its confinement in a complex magnetic topology. They calculated the footprints of evolving quasi-separatrix layers (QSLs 146 ) during the flare reconnection, which are found to coincide rather well with the structure and evolution of the observed flare ribbons. By careful analysis of different components of the Lorentz force (Figure 8), they revealed for the first time that a particular component, resulting from the MFR's radial magnetic field, which is defined in the local cylindrical coordinates of the MFR and is related to the non-axisymmetry of the MFR's cross-section, essentially constrains the erupting MFR. This differs from the prediction of TI, in which the major confining force should come from the external strapping field, and such new results demonstrate that a data-driven model can disclose the complexity of physics in a realistic eruption that has not been considered in traditional theories. The data-driven model has also been shown to be able to simulate successive eruptions from the same region. 38

SUMMARY AND OUTLOOK
To summarize, major developments of data-driven modeling of coronal magnetic fields have been witnessed over the last 10 years. A variety of models have been proposed, with different assumptions and implementations of the controlling equations, the numerical solver, as well as the initial and boundary conditions to drive the models. These models have been applied to simulation of the long-term evolution and fast eruptions of complex coronal magnetic fields with various topologies. The results have well demonstrated that data-driven simulations are able to achieve an unprecedented high degree of realism exceeding non-data-driven simulations, in particular by reproducing the eruption process with a complexity that is beyond the scope of idealized or theoretical models. This has been shown by the striking resemblance of the magnetic topology in the simulated eruptions with the observed erupting coronal loops and filaments as well as the complex morphology and evolution of the flare ribbons. Moreover, by recovering a continual evolution that allows a "memory" of magnetic topology and build-up of electric currents and free magnetic energy, data-driven models maintain a high degree of self-consistence compared with other models. Data-driven modeling helps us to find the key mechanism that initiates the eruptions and reveal more in depth the physical behaviors of eruption that challenge theoretical and idealized models, which has greatly extended our understanding of the coronal magnetic dynamics. Although different eruption scenarios can be tested with simulations of other types, only data-driven simulations can be used to identify the actual scenario in real events. For example, many NLFFF analyses 84,157-160 of coronal fields show pre-flare MFR, suggesting that the MHD instabilities of the MFR could trigger the eruption (by studying the twist number and decay index), but only data-driven MHD simulations can determine whether these instabilities exist and their role in producing eruptions. Some unexpected results have already been found. For instance, a Lorentz force component resulting from an MFR's radial magnetic field could result in failure of the MFR to erupt, which is not considered in the theory of TI. Another finding is that, when reconnection begins, the eruption is mainly driven by the ongoing reconnection rather than TI, which explains why eruptions can still succeed in a background field with a saddle-like profile of the decay index.
Next we provide a summary of future developments and applications of these models.

Joint assessment of different codes
As in the case of NLFFF extrapolation, for which a variety of codes exist and joint testing of these codes is important to guide code improvements, 147,149,161 different data-driven models should also be compared jointly; for example, using some ground-truth datasets (but of course from other types of simulations). Such a comparative study can help to identify the key factors that make performance of the various models different. Toriumi et al. 162 first examined four different datadriven models 23,35,43,44 using data from an idealized flux emergence simulation (FES) of a twisted magnetic flux tube emerging from the convection zone into the corona. The photospheric magnetic and velocity fields produced by the FES are input to drive the different data-driven models as part of their bottom boundaries. It was found that, although all of these data-driven models successfully reproduced a flux rope structure, the quantitative discrepancies are considerable. Even though some of the models can reproduce the magnetic field energy and relative helicity in the corona quite well, the 3D magnetic configuration is still rather different from the ground-truth data (see Figures 4 and 5 in Toriumi et al. 162 ). The discrepancies among the models were attributed mainly to the highly non-forcefree input photospheric field and the modeling constraints, such as the treatment of background atmosphere, the bottom boundary setting, etc.
Indeed, in a following work, Jiang and Toriumi 163 tested the DARE-MHD model with the bottom boundary driven by three sets of data from the same FES sliced at incremental heights, which correspond to the photosphere, the chromosphere, and the base of the corona, respectively. They first compared the relative magnitudes of the Lorentz force contained in these datasets using the ratio of the integrated Lorentz force (and torque) to the integrated magnetic pressure gradient force (and its torque). This ratio, or normalized force (and torque), for a perfectly force-free magnetogram is 0 (a criterion of below 0.1 is commonly accepted to be force free [164][165][166], whereas, for an extremely non-force-free field, it is close to 1, meaning that the magnetic pressure force and the tension force are so unbalanced that the net Lorentz force is comparable with one of its components, the total magnetic pressure force. As found by Jiang and Toriumi 163 , for the photospheric data of the FES, the ratio (for force and torque) is close to 1 in the whole time evolution, meaning that the photosphere data have a very strong Lorentz force. For the data sliced at two higher levels, the normalized Lorentz force and torque in the early phase of emergence are still large (close to 1), but decrease substantially (to around 0.2) with emergence of the magnetic flux, indicating that, at these levels, the fields are relaxing to a more force-free state. As a result, the DARE-MHD using the photosphere data fails to reproduce the coronal magnetic field, whereas at the two higher levels above the photosphere, the DARE-MHD yields results in much better agreement with the FES. This confirms that the Lorentz force in the boundary data is a key issue affecting the results of the DARE-MHD model, and the key reason is that the bottom surface in the DARE-MHD model (and in many other models) is assumed to be the coronal base with a plasma number density of around 10 9 cm -3 , which is many orders of magnitude smaller than that of the photosphere (about 10 17 cm À3 ). Thus, if the strongly forced photospheric magnetic field is entered into a model that uses typical The gray line represents the net force. The orange, purple, blue, red, and cyan lines represent the different components, including the hoop force (F H ), non-axisymmetry induced forces (F N1 and F N2 ), strapping force (F S ), and tension force (F T ). The black vertical line indicates the location of the axis of the MFR. F N1 changes its direction during the eruption process. It is initially upward but becomes downward in the end, constraining the MFR eruption. Image reproduced from Zhong et al. 39 settings for the atmosphere in the corona, such a strong Lorentz force cannot be balanced by the tenuous plasma. It can induce largely spurious plasma motions, which, in turn, amplify the magnetic field by the magnetic induction equation and could make the evolution run away. 162,163 The other two data-driven MHD models, developed by Guo et al. 35 and Hayashi et al. 43 and examined by Toriumi et al., 162 use much denser plasma (with a number density of 10 13 cm À3 and 10 17 cm À3 , respectively) at the lower boundary and can avoid such an inconsistency (by keeping the induced motions at a reasonable level of speed). This emphasizes the need for incorporating a more realistic atmospheric stratification from the photosphere, chromosphere, and transition region to the corona (if the photosphere magnetic field is the only data that can be used). With this, the data-driven model can also be used to probe the physical mechanisms behind the reverse feedback of coronal eruption to the photosphere; for example, the rapid, permanent change of photospheric magnetic fields as often observed through flares 167,168 and the sudden reversal of rotation in a sunspot during a solar flare. 169,170 On the other hand, the observed photospheric field for flux-emerging ARs from SDO/HMI has been found to be much closer to force free than the FES data, as demonstrated recently by Duan et al. 171 They conducted a statistical study of the normalized Lorentz forces and torques with a substantially large sample of 3,536 vector magnetograms from 51 ARs taken by SDO/HMI over the period of 2010-2019. These ARs were selected as having significant flux emergence during their passage on the solar disk within ± 45 in longitude from the Sun's central meridian. The normalized Lorentz forces and torques for all vector magnetograms are, on average, very close to 0.1, which is much smaller than that from the simplified FES. Thus, on one hand, in future investigations, one should consider using realistic convective FES models that could yield more relaxed emergence 172,173 as the ground-truth data to perform a more insightful assessment of the different data-driven models. Or, in comparison with other data-driven models, the models that include the atmosphere of only the corona should be provided with magnetic field data in the higher and more force-free solar atmosphere, especially considering that some of these models adopt an NLFFF as the initial condition. On the other hand, such an inconsistency, as found by the FES and the observed magnetogram (i.e., the latter is much more closer to force free than the former) calls for more realistic flux emergence models 174 and even solar dynamo models 175 with constraint from the observed photospheric magnetic field.

Coupling of local and global models
Currently, most of the data-driven models for eruptions are implemented in a local Cartesian box, which, unavoidably, has numerical boundaries at the side and top of the computational volume, and it has long been a big problem to specify the conditions on these boundaries so that they can be numerically stable and non-reflecting. To avoid this difficulty, many modelers choose to simply fix the variables on the numerical boundaries and stop their simulation before the disturbance induced by eruption reaches these boundaries. 35,39,109 There are also models using periodic conditions for the side boundaries, which, however, may result in unreasonably distorted structures of the simulated magnetic field and plasma. 176 Future data-driven simulations of eruptions need to be extended to the global corona and incorporate solar wind by utilizing full-disk vector magnetograms (rather than only LOS synoptic maps) to study how an initiated CME evolves and interacts with the large-scale structures in the local corona and in the interplanetary space. This, on one hand, can avoid the sub-Alfvénic numerical boundaries, and more importantly, on the other hand, can provide a more realistic evolution of CME in the heliosphere than most current simulations of CME propagations that commonly employ an artificially inserted, force-unbalanced MFR in the corona to generate an eruption [177][178][179][180][181] and will definitely be a big step forward for developing sophisticated model of space weather, which is driven by solar eruptions. By coupling the local and global models, it is also important to examine whether the eruptive structure (e.g., flux rope) formed in the local frame is confined or has access to outer space. For instance, DeRosa and Barnes 182 surveyed the large-scale magnetic flux overlying the ARs of a number of X-class flares (including eruptive and confined ones) using the potential field source surface (PFSS) model. 15 They found that, in general, the ARs of eruptive events are located closer to the open flux than those of confined ones. These kinds of studies using datadriven models may provide important insights into how the eruptive structures interact with the different background fields and result in successful or failed eruptions. A few efforts have already been made to develop such coupled models. For example, Wu et al. 183 entered the eruptive magnetic structure of AR NOAA 11283 obtained from Jiang et al. 101 into a global coronal and solar wind model to track the propagation of the eruption. Hayashi et al. 184 used the output from the magnetofrictional model (spherical version) of Cheung and DeRosa 23 to specify the inner boundary conditions of their global heliosphere model and simulated an eruption from AR NOAA 11158. More quantitative comparisons between the simulation results with coronagraph observations and in-situ data within the corona are required; for example, from the Parker Solar Probe 185,186 and Solar Orbiter, 187 to validate and constrain the data-driven models.
At the other end of the spectrum, data-driven MHD simulations for eruptions should also be coupled with very small-scale photospheric magnetic field evolutions, as observed with very high spatial resolution but in a small field of view. This is because there are small-scale processes playing key roles in producing eruptions. For example, some high-resolution observations 188,189 show that the emergence of very-small-scale magnetic flux reversed to the large-scale field of a pre-existing sheared arcade, the so-called reversed magnetic shear mechanism, 139,188 could lead to a global eruption, and this mechanism has yet to be investigated by data-driven modeling. Another example is magnetic flux cancellation, in which small-scale magnetic flux elements of opposite polarities meet at the PIL and cancel each other; this is believed to be a key mechanism of building up pre-eruption coronal MFRs. Because the nature of flux cancellation is still elusive, this needs to be carefully studied with data-driven simulations driven by high-resolution data focusing on the small-scale magnetic structures near the PIL. High-spatial-resolution vector magnetograms obtained by, e.g., the Daniel K. Inouye Solar Telescope (DKIST 190 ), could provide boundary data for such small-scale MHD simulations. High-temporal-resolution data are also indispensable for simulating fast flux emergence events. 72 In addition to the spatial resolution and temporal cadence, vector magnetic fields also have to be accurate enough to maintain the credibility of data-driven simulations. Vector magnetic field observations often suffer from two major problems, 180 ambiguity and the relatively low accuracy of the transverse component compared with the LOS component. Significant progress to resolve these problems was made by launching the Solar Orbiter. 187 The combination of vector magnetic field observations from two perspectives, such as SDO and Solar Orbiter, could resolve the 180 ambiguity 191 and increase the precision of their transverse components.

Including more physics
Most data-driven models as discussed are focused on magnetic fields, whereas the plasma is not dealt with realistically (either for initial conditions or in model equations), especially in the magnetofrictional and the zero-b models. By discarding the thermal pressure, these simplified models might fail when there is a fast reconnection in the field, where thermal pressure could play an important role in the dynamics in the weak-field region of magnetic field dissipation. Further improvements should be devoted more to the full MHD and even the radiative MHD models. On the one hand, more physics need to be included in the energy equation by considering the effects of coronal heating, thermal conduction, and radiative cooling. On the other hand, more realistic initial and boundary conditions are required to be used for not only the magnetic field but also the plasma density and temperature, which, for example, could be determined using differential emission measure tomography. 192 With these advancements, the models can obtain the capability to simulate the wealth of emission structures (such as EUV coronal loops, filaments, and flare emissions) of the coronal field that can be directly compared with observations, which, in turn, gives more constrains to the models.

Achieving higher accuracy of computation
The value of numerical resistivity is also a problem for modeling magnetic reconnection. Currently, most data-driven models have numerical resistivity with values much higher by orders of magnitude than the realistic coronal value, and, thus, slow reconnection is often triggered too easily; in other words, the current layer cannot be sufficiently thin to let a fast reconnection set in. For example, a recent simulation 97 shows that the magnitude of resistivity can substantially affect the result of eruption-initiation simulations, which is associated with a long-standing question regarding the initiation of eruptions. It has been well established by many MHD simulations [193][194][195] in two dimensions (or translationally invariant geometries) with high accuracy that a magnetic arcade sheared continuously will asymptotically approach an open field 196,197 containing a thin current sheet, and the system experiences a global eruption when reconnection sets in at the current sheet. However, this simple and efficient scenario of eruption initiation has never been reproduced in fully 3D MHD simulations with sufficiently high accuracy until recently by Jiang et al., 97 and the key reason why previous 3D simulations failed is that their models have an unreasonably large resistivity (in the line-tied bottom boundary, the computational volume, or both) so that the thin current sheet cannot form (because the current layer dissipates quickly by the large resistivity). Therefore, future data-driven models should consider improving numerical accuracy and resolution to reach a sufficiently small resistivity and, eventually, to bridge the gap between the numerical Lundquist number of only a few 10 2 $ 10 3 to a realistic one of around 10 13 .

Improving computational efficiency
A major limitation for application of data-driven comprehensive models that should include all aforementioned improvements is the large amount of computing resources required. According to the authors' experience, to simulate the evolution of a typical size AR (say, 300 Mm in all three directions) for a typical timescale of a few days with a grid resolution that matches the HMI data (720 km), an MHD model with isothermal simplification needs months of computing time even when parallelized with a medium number of CPUs (for example, 100 processors with a frequency of 3 GHz). To reduce the computing time, some modelers 44,151,152 chose to speed up the cadence of feeding the observed data into the model by tens of times. This is reasonable by taking advantage of the fact that the photospheric evolution speed (with a typical velocity of 0.1-1 km s À1 ) is slower than the coronal evolution speed (of a few Mm s À1 ) by about three orders of magnitude. Speeding up the driving process is also necessary, considering that the numerical resistivity is much larger than the realistic one, and, thus, the global diffusion speed of the free magnetic energy could be faster than the energy injection rate from the bottom boundary when using the realistic timing of driving. In this case, no free energy (and current) can be accumulated in the corona, and no eruption can be initiated. Although speeding up photospheric driving can be acceptable in the non-eruptive time duration, the coronal magnetic field may still not relax sufficiently to reproduce the quasi-static evolution in the real corona, while during eruptions, it may not be consistent with the sped-up evolution of the photospheric field. Therefore, solving all these problems demands more accurate numerical solving of the controlling equations and faster computing techniques developed for hardware (for example, using GPU computing) and software. The next generation of data-driven modeling should be efficient enough to match the fast-developing observation capability and eventually have computing power faster than real time and an accuracy high enough for prediction of solar eruptions and for use as a key module for space weather forecasting.