Simulating the Coronal Evolution of Bipolar Active Regions to Investigate the Formation of Flux Ropes

The coronal magnetic field evolution of 20 bipolar active regions (ARs) is simulated from their emergence to decay using the time-dependent nonlinear force-free field method of Mackay et al. A time sequence of cleaned photospheric line-of-sight magnetograms, that covers the entire evolution of each AR, is used to drive the simulation. A comparison of the simulated coronal magnetic field with the 171 and 193 A observations obtained by the Solar Dynamics Observatory (SDO)/ Atmospheric Imaging Assembly (AIA), is made for each AR by manual inspection. The results show that it is possible to reproduce the evolution of the main coronal features such as small- and large-scale coronal loops, filaments and sheared structures for 80% of the ARs. Varying the boundary and initial conditions, along with the addition of physical effects such as Ohmic diffusion, hyperdiffusion and a horizontal magnetic field injection at the photosphere, improves the match between the observations and simulated coronal evolution by 20%. The simulations were able to reproduce the build-up to eruption for 50% of the observed eruptions associated with the ARs. The mean unsigned time difference between the eruptions occurring in the observations compared to the time of eruption onset in the simulations was found to be ~5 hrs. The simulations were particularly successful in capturing the build-up to eruption for all four eruptions that originated from the internal polarity inversion line of the ARs. The technique was less successful in reproducing the onset of eruptions that originated from the periphery of ARs and large-scale coronal structures. For these cases global, rather than local, nonlinear force-free field models must be used. While the technique has shown some success, eruptions that occur in quick succession are difficult to reproduce by this method and future iterations of the model need to address this.


Introduction
The solar corona is highly complex in nature. The source of its complexity is largely due to the presence of magnetic fields that are generated in the tachocline (Spiegel and Zahn, 1992): a region close to the base of the convection zone (Charbonneau, 2010(Charbonneau, , 2014. When magnetic flux tubes at the base of the convection zone become unstable to buoyancy (Parker, 1955;Zwaan, 1985) they rise and the magnetic field breaks through the solar surface manifesting itself as an active region (AR) in the photosphere. The magnetic flux emerges in a non-potential state (Leka et al., 1996) and is further modified by the action of photospheric flows. This results in free magnetic energy being available to drive solar eruptive phenomena.
ARs are the source of a wide range of atmospheric solar activity and the type and level of activity is dependent on the evolutionary stage of the AR (for a review on AR evolution see van Driel-Gesztelyi and Green 2015). As a result, it is important to understand the structure and evolution of the magnetic field of an AR over its entire lifetime, from emergence to decay.
It is currently difficult to measure the magnetic field in the corona and extreme ultraviolet (EUV) observations of AR coronal loops can only provide indirect and limited information of the coronal structure of ARs. An alternative approach, for the analysis of the coronal structure of ARs, is to construct a model of the coronal magnetic field by using the photospheric magnetic field as the lower boundary condition. This approach relies on the approximation that the corona, a low plasma-β environment that mostly remains in equilibrium, is "force-free". This means that the coronal magnetic field must satisfy the criterion of j × B = 0 where j = αB. In the case of nonlinear force-free (NLFF) fields the torsion parameter α is a scalar function that can vary as a function of position, but must remain constant along magnetic field lines.
There are numerous NLFF field techniques that can be used to generate models of the coronal magnetic field. These NLFF field models can be divided into two categories: models that are static or time-dependent. Static models either use a vector magnetogram as the lower boundary condition and extrapolate the NLFF fields into the corona (e.g. Schrijver et al. 2006;De Rosa et al. 2009;Canou and Amari 2010;Wiegelmann and Sakurai 2012;Jiang et al. 2014), or they take an initial coronal field, which is either a potential or linear force-free (LFF), and evolve this field into a NLFF state. The latter approach can make use of the magnetofrictional relaxation technique (Yang, Sturrock and Antiochos, 1986) to generate a static model of the magnetic field of an AR. Examples of static modelling using magnetofrictional relaxation include the magnetofrictional extrapolation method of Valori, Kliem and Keppens (2005) and the flux rope insertion method (van Ballegooijen, 2004;Bobra, van Ballegooijen and DeLuca, 2008;Savcheva et al., 2012;Yardley et al., 2019). The extrapolation methods mentioned above produce a coronal field model at a single snapshot in time.
A series of independent, static extrapolations may be produced but there is no direct evolution from one extrapolation to the next.
The magnetofrictional relaxation technique can also be used as a simulation method to construct a continuous time-dependent series of NLFF fields. In this case, the normal component of the magnetic field is specified along with an initial field and a time series of horizontal boundary motions. The resulting coronal structures are due to the applied boundary motions injecting non-potentiality into the corona over timescales of hours or days. The coronal field, which is in non-equilibrium is then relaxed back to a NLFF field equilibrium using magnetofrictional relaxation. This has been applied to global simulations (Mackay and van Ballegooijen, 2006a,b) where a flux transport model is applied at the photospheric boundary or to simulate AR evolution using a time series of lineof-sight (LoS) magetograms (Mackay, Green and van Ballegooijen, 2011;Gibb et al., 2014) or more recently vector magnetograms (e.g. Pomoell, Lumme and Kilpua 2019).
In the recent study by Yardley, Mackay and Green (2018b) a continuous time-dependent series of NLFF field models of AR 11437 were created using the time-dependent NLFF field method of Mackay, Green and van Ballegooijen (2011). Photospheric LoS magnetograms from the SDO/Helioseismic Magnetic Imager (HMI) instrument were used as lower boundary conditions to drive the simulation and continuously evolve the coronal field through a series of NLFF equilibria. When the results from the simulation were compared to SDO/AIA observations it was found that the simulation was able to capture the majority of the characteristics of the coronal field evolution. Flux ropes that formed in the simulation showed signatures of eruption onset for two out of three of the observed eruptions, approximately 1 and 10 hrs before the eruptions occurred in the observations. A parameter study was also conducted to test whether varying the initial condition and boundary conditions along with the inclusion of Ohmic diffusion, hyperdiffusion, and an additional horizontal magnetic field injection at the photosphere affect the coronal evolution and timings of the eruption onset. The results showed that the coronal evolution and timings of eruption onset were not significantly changed by these variations and inclusions, indicating that the main element in replicating the coronal field evolution is the Poynting flux from the boundary evolution of the LoS magnetograms. AR 11437 is also included in this current study.
In this paper, we extend the set of simulations carried out in Yardley, Mackay and Green (2018b) of a single AR by simulating the coronal magnetic field evolution of 20 bipolar ARs. The observational analysis of the same set of bipolar ARs was conducted by Yardley et al. (2018a) in order to probe the role of flux cancellation as an eruption trigger mechanism. The study of Yardley et al. (2018a) analysed both photospheric and coronal observations taken by SDO over the entire lifetime of the ARs. Through simulating a much larger sample of ARs we can obtain more general results than those found in Yardley, Mackay and Green (2018b), which only considered a single region (AR 11437). We aim to determine whether the simulation of a series of NLFF fields using the magnetofrictional technique can capture the coronal evolution and also the build-up phase that brings the coronal field to the point of eruption. The analysis carried out here is similar to that of Yardley, Mackay and Green (2018b) in which the NLFF field method was tested. However, due to the large-scale analysis of 20 ARs the results are presented in less detail than those given in Yardley, Mackay and Green (2018b).
The outline of the paper is as follows. Section 2 outlines the observations including the criteria for AR selection, coronal evolution and eruptions produced by each AR. Section 3 describes the technique used to simulate the coronal field including the lower boundary conditions used. Results from the simulations can be found in Section 4, which includes simulations using the simplest initial and boundary conditions and also with the inclusion of additional effects. Section 5 discusses the results and Section 6 provides a conclusion to the study.

AR Selection
The 20 ARs presented in Yardley et al. (2018a) are the same regions used in this study. We now briefly summarise the data selection method used by Yardley et al. (2018a) to identify and select these ARs and refer the reader to that paper for more details on each region. ARs were selected using the following criteria: i) The ARs must be bipolar and have low complexity. The regions must have two dominant photospheric magnetic polarities with no major mixing of the opposite polarities. ii) The ARs must be isolated with minimal interaction occurring between the AR and other ARs or the background quiet Sun magnetic field. iii) The ARs must be observable from their first emergence and form east of central meridian. This allows the full evolution from emergence to decay to be simulated during disk transit. iv) The ARs first emergence must be no more than 60 • from central meridian as instrumental effects become increasingly significant at large centre-to-limb angles.
These selection criteria led to a sample of 20 ARs being chosen during the HMI era, spanning a time period from March 2012 to November 2015. All ARs, apart from AR 11867, were monitored during their flux emergence and decay phases, which included dispersal and flux cancellation. AR 11867 remained in its emergence phase during the time period studied and did not exhibit flux cancellation at its internal PIL.
Representative AR examples are given in Figure 1 with Supplementary Movie 1 showing the full evolution of AR 11446. Table 1 provides summary information of AR locations, photospheric flux evolution, and observed eruption times taken from Yardley et al. (2018a). Photospheric flux values were obtained using the 720 s data series (Couvidat et al., 2016) generated by the Helioseismic Magnetic Imager (HMI) (Schou et al., 2012) (ARs 11437, 11446 & 11680). The images show each AR at the time of the peak unsigned magnetic flux measurement, where unsigned refers to half the total absolute positive and negative flux. The saturation levels of the images are ± 100 G with white (black) representing positive (negative) photospheric magnetic field. As an example, the entire photospheric field evolution of AR 11437 can be seen online in Supplementary Movie 1.

Coronal Evolution and Eruptive Activity
The observed coronal evolution of each AR was analysed in Yardley et al. (2018a) in order to identify the time and location of any eruptions. These ejections are referred to as eruptions as opposed to CMEs because the coronal signatures in the EUV data are relatively subtle and most do not show any clear evidence of a CME in the white-light coronagraph data. This implies that they are either confined/failed eruptions or are ejective but have a low plasma density. The coronal evolution was monitored using both 171 and 193Å images taken by the Atmospheric Imaging Assembly (AIA; Lemen et al. 2012) on board SDO. AIA provides full-disk observations with a high spatial and temporal resolution of 1.5" and 12 s, respectively. At least two or more of the following coronal signatures were used to identify the occurrence of an eruption: i) the eruption of a filament or an EUV loop system, ii) the rapid disappearance of coronal loops and post-eruption arcade formation (flare arcade), iii) flares and flare ribbons, iv) and/or coronal dimmings.
As detailed in Yardley et al. (2018a) the eruptions were then categorized into the following types to investigate which eruptive structures might have formed as a consequence of flux cancellation: i) Internal PIL events are the eruption of a low altitude structure originating along the internal PIL of the AR. ii) External PIL events are the eruption of a low altitude structure originating along an external PIL that is formed between the periphery of the AR and the magnetic field of the quiet Sun. iii) High altitude events are the eruption of a high altitude structure which cannot be associated with an internal/ external PIL (which are at low altitude).
In total, 24 eruptions were observed, with 13 of the 20 ARs producing at least one ejection. Eight of these ARs produced low corona events originating from either the internal or external PIL and the other five produced high altitude events. Two of the eruptions were observed as a CME in the LASCO/C2 coronagraph data. There were also four B/C GOES class flares associated with four ARs that did not occur at the time of the eruptions.
For examples of the different event categories see Figure 1 in Yardley et al. (2018a). The timings of these events, which are also taken from Yardley et al. (2018a), are given in Table 1.

Coronal Magnetic Field Evolution
The NLFF field method of Mackay, Green and van Ballegooijen (2011) is applied to SDO/HMI LoS magnetograms to simulate the evolution of the coronal magnetic field of each AR. A key element of this method is that the magnetic field evolves through a continuous time series both at the photosphere and in the coronal volume where flux is preserved. Therefore, the coronal magnetic field evolution can be analysed. When using our method we do not apply any additional observational constraints such as the use of EUV coronal images, rather the solution obtained at any one time is purely based on the initial field, the applied boundary motions and any additional coronal physics (see Section 4.2).
This technique has been previously tested on AR 11437 (Yardley, Mackay and Green, 2018b), one of the ARs also included in this study. Therefore, the quantitative analysis that has previously been carried out for AR 11437 will not be described in this paper. Here, we present the overarching results from the qualitative analysis of 20 bipolar ARs, where each AR has been studied using the methodology described in Yardley, Mackay and Green (2018b).
A time series of NLFF fields is generated using HMI LoS magnetograms for each lower boundary condition (see Section 3.2). The HMI LoS magnetograms are cleaned and re-scaled before the simulations are carried out. The clean-up procedure includes time-averaging, low magnetic flux value removal, removal of small-scale magnetic elements, and if required, flux balancing. This procedure ensures that the large-scale AR evolution is kept but small-scale quiet Sun elements and random noise are removed (see Appendix A for more details).
In the simulation, the evolution of the 3D magnetic field B is described by where A represents the magnetic vector potential, B = ∇ × A is the magnetic field, and v is the magnetofrictional velocity. The magnetofrictional relaxation technique of Yang, Sturrock and Antiochos (1986) is employed to ensure that the coronal field is evolved through a series of force-free equilibria. Therefore, the magnetofrictional velocity inside the computational box takes the form where ν is the friction coefficient and j = ∇ × B. The coefficient of friction ensures that as the magnetic field is perturbed by motions at the boundary, the field remains close to a force-free equilibrium in the corona. A cartesian staggered grid is used to carry out the computations to obtain second-order accuracy for A, B, and j.

Photospheric Boundary Conditions
To be able to simulate the full evolution of the bipolar ARs we use the full disk HMI 720s LoS magnetograms (hmi.M 720s series). For each AR, we use a time sequence of LoS magnetograms with a chosen cadence of 96 minutes. We create cut-outs of the magnetograms centred on each AR and apply clean-up processes to the time series of partial disk magnetograms (see Appendix A). We use LoS magnetograms in this study as we want to simulate the full evolution of ARs from emergence to decay. We would also like to quantify how well this computationally efficient modelling technique that uses the LoS magnetograms performs in simulating the coronal evolution of a large number of ARs. Regarding the medium cadence used, prior to the present study, we have conducted a number of investigations varying the cadence of the HMI magnetograms from 12 minutes to 3 hours (Gibb, 2015) and have found very similar results. Therefore, we have chosen to use a medium cadence of 96 minutes as it is sufficient to capture the large-scale evolution of the ARs. Also, any future L5 space weather mission is likely to have a cadence more comparable to that of the Michelson Doppler Imager (MDI) rather than the cadence presently provided by HMI.
Initially, each simulation is run using a relatively simple set-up. That is, a potential magnetic field is used as the initial condition along with either a closed or open boundary at the top of the computational volume. The simulation results are then compared with the observations to determine whether there is a good agreement between the two. This is assessed by comparing the evolution of the simulated coronal field to the coronal evolution in SDO/AIA 171 and 193Å observations by visual inspection and using qualitative scoring criteria given in Section 4.1. If the simulation results do not provide a good fit to the observations then the simulation is re-run varying a number of terms one-by-one. First a LFF field initial condition is used then a variety of additional physical effects are included in succession until a better fit is achieved (see Section 4.2 and also the method of Yardley, Mackay and Green 2018b). For the present simulations we only use a potential or a LFF field as the initial condition for the simulations. The ARs modelled in this study are young ARs, with the majority emerging at a centre-to-limb angle around 60 • longitude. Due to the large distance from central meridian the vector magnetograms (where they may exist) contain significant errors where these errors could introduce spurious results in the simulation. Therefore, using an initial NLFF field condition is currently beyond the scope of this paper but this will be considered in a future study.
The simulations use the cleaned LoS magnetograms (see Appendix A), which have a been scaled to a lower resolution of 256 2 , as the lower boundary conditions. The original size of the magnetograms depends upon the size of the AR but the LoS magnetograms are always larger than 256 2 . To take into account boundary effects, the magnetograms are also re-scaled to fill 60-70% of the area of the bottom of the computational box. The simulation generates a continuous series of lower boundary conditions using the corrected LoS magnetograms that are designed to replicate, pixel by pixel the LoS magnetograms, every 96 minutes.
The series of cleaned magnetograms give the prescribed distribution of B z on the base. Hence the horizontal components of the vector potential (A xb , A yb ) are determined on the base for each discrete time interval of 96 minutes by solving for the scalar potential φ, where A = ∇ × (φẑ). To specify the evolution of B z on the base in terms of A xb and A yb between the prescribed distributions the rate of change of the horizontal components of the magnetic vector potential and therefore an electric field is determined. To evolve A xb (t) and A yb (t) to A xb (t+1) and A yb (t + 1) we assume that the process is linearly applied between each discrete time interval t and t + 1, where t represents the discrete 96 minute time index. Therefore, the horizontal components (A xb , A yb ) are linearly interpolated between each 96 minute time interval to produce a time sequence that is continuous between the observed distributions. Thus, every 96 minutes the simulated photospheric field identically matches that found in the cleaned observations. By using this technique, we are effectively evolving the magnetic field from one fixed magnetogram to the next. Also, undesirable effects such as the pile-up of magnetic flux at sites of flux cancellation and numerical overshoot do not occur. As the surface field evolves in this manner it injects electric currents and free energy into the coronal field, which responds through Equation 1. By using this numerical method it means that there are two timescales involved in the lower boundary condition evolution. The first timescale is due to the 96 minute time cadence of the observations and the second is the linear evolution timescale. The second timescale is introduced to advect the photospheric magnetic polarities between the observed states, inject Poynting flux into the corona and to relax the coronal field. The method applied to interpolate the boundary magnetic field is very similar to Gibb et al. (2014) and Yardley, Mackay and Green (2018b) however, to satisfy the Courant-Friedrichs-Levy (CFL) condition the timestep is determined from the minimum cell crossing time for the magnetofrictional velocity or the diffusion terms and its maximum is equal to a fifth of this value.
Within the simulations the initial condition satisfies the Coulomb gauge. In addition to this, during the evolution of the field between the fixed points given by the magnetograms, we also maintain the Coulomb gauge. This is carried out numerically by including a ∇ · A term, which does not effect the value of the magnetic field in the simulations. The complete description of this process can be found in Mackay and van Ballegooijen (2009), Mackay, Green and van Ballegooijen (2011) and references therein.

Magnetic Field Evolution
The simulated coronal field evolution of the 20 bipolar ARs will now be discussed for the simplest case where a potential field is used as the initial condition and the top boundary of the computational box is closed. To determine whether the simulated coronal evolution is able to capture that of the real Sun in each AR, the simulated field is compared to the SDO/AIA 171 and 193Å plasma emission structures by manual inspection. The main coronal features that are used to make the comparison between the observed coronal structure and simulated coronal magnetic field of each AR include small-and large-scale coronal loops, filaments and sheared structures. The 171 and 193Å wavebands are used for the comparison as the evolution of coronal loops, filaments and sheared structures are well captured in these wavebands compared to the other AIA wavebands. These wavebands are also the primary wavebands that were analysed in the observational study of Yardley et al. (2018a). The simulated magnetic field and observed coronal plasma emission structures are then compared at various times (roughly once per day, see Figure 2) throughout the evolution of each AR.
The simulation results are also analysed to determine whether or not there is a good agreement between the timings and location of the ejections seen in the observations and the corresponding signatures of eruption onset in the simulations. The following criteria are used to assign a score to quantitatively describe the level of agreement between the simulations and observations: i) Score 1: If the simulation is able to reproduce the main coronal features (smalland large-scale loops, filaments and sheared structures) for the majority of the AR evolution then there is deemed to be a good match between the observations and simulations. If an eruption is observed to originate from the AR, then the simulation must be able to successfully model the build-up to the eruption within a ±12 hr time window pre-or post-observed eruption time. If there are multiple observed eruptions then the simulation must be able to successfully follow the build-up to eruption for the majority of the eruptions associated with the AR.
ii) Score 2: Some of the coronal features (small-and large-scale loops, filaments and sheared structures) that are seen in the observations are reproduced by the simulation for most of the AR evolution. Therefore, the match between the coronal features present in the observations and the simulations is deemed to be acceptable. If one or multiple eruptions are observed to originate from the AR, the build-up phase may or may not be followed by the simulation for any eruption. iii) Score 3: A minority or none of the coronal features (small-and large-scale loops, filaments and sheared structures) seen in the observations are reproduced for most of the AR evolution. Therefore, the evolution of the simulated coronal field is deemed not to match the observed coronal evolution. The simulation fails to model the build-up to eruption for any observed eruptions associated with the AR.
An example AR for each of the scoring criteria is shown in Figure 2, which compares the observed coronal evolution (odd rows) to the simulated coronal evolution (even rows). The first example shows AR 11437 (Score 1), where the sheared J-shaped structure, small-and large-scale coronal loops that are present in the observations are captured by the simulation for the majority of the AR evolution (see black arrows in Figure 2). The simulation is also able to replicate the build-up to the point of eruption for 2 out of 3 of the observed eruptions. The signatures of eruption onset in the simulations are discussed in the next paragraph. The second example shows AR 12455 (Score 2), where the simulation is able to reproduce the structure of the small-and large-scale coronal loops although, the match to the observations is better in the northern part of the AR compared to the south (black arrows in Figure 2). There are no eruptions observed to be associated with this AR. Finally, for AR 12229 (Score 3) the simulation is unable to produce the structure of the small-and large-scale loops seen in the observations of the AR. The eruption onset signatures, which indicate that a loss of equilibrium in the simulation has occurred, are not present for any of the four eruptions observed to originate from this AR.
The simulations carried out, focus on modelling the build-up of non-potential magnetic fields and flux ropes within ARs. We do not try to reproduce and follow the dynamics of the observed eruptions as full magnetohydrodynamic (MHD) simulations are required to do this (e.g. see Rodkin et al. 2017). Therefore, to determine whether the simulations successfully follow the build-up to eruption, the simulated coronal field evolution was examined for signatures of eruption onset. The signatures present in the simulations that indicate the build-up to an eruption include: i) a flux rope rising, which subsequently reaches the top or side boundaries of the computational box indicating that a loss of equilibrium has occurred. ii) Reconnection occurring underneath the flux rope which leads to small, more potential loops forming beneath the flux rope similar to the post-eruption (flare) arcades that are visible in the observations.
These signatures of eruption onset in the simulations must occur at the same location and timings as those identified in the observations. The simulation results are analyzed in a time window of ≈12 hrs pre-and post-observed eruption for the above signatures of eruption onset. The signatures of eruption onset in the simulation of AR 11437 are shown in Figure 3. In this case, a flux rope, which has formed along the internal PIL, rises in the domain and reconnection occurs underneath the flux rope. This leads to small, more potential loops forming below the flux rope axis. Eventually the flux rope reaches the side boundary of the domain. A similar scenario is seen in the observations where a sheared structure and post-eruption loops that form underneath this structure are observed at the same location as in the simulations. For the simplest case, where the coronal evolution of each of the 20 bipolar ARs is simulated using a potential field initial condition and a closed top bound- ary, the results (see Table 2) are as follows. The NLFF field method is able to capture the majority of the coronal structure for ten ARs, a reasonable amount of the structure for six ARs, and little or no structure for four ARs (see Table 3). Therefore, the method is able to capture a reasonable amount of the structure for 80% of the AR sample, and failed to capture the structure for 20% of the ARs.
In total, the simulations are able to successfully follow the build-up to eruption in a ≈12 hr time window prior to or post-eruption for 12 out of the 24 observed eruptions. The time difference between eruption onset in the simulations compared to the time determined from observations is given in Figure 4 for each AR. The time of eruption onset in the simulation is determined by using the time halfway between the time step where the signatures of eruption onset in the simulation have been identified, and the previous time step where there are no signatures of eruption onset. By time step we are referring to the primary timescale of the simulation that is set by the cadence of the magnetograms, which in this case is 96 minutes. The time of the eruption onset identified from the simulation is then compared to the eruption time taken from the observations to give the time difference. The mean time difference between the initiation of the eruption in the simulations compared to the observations is ≈5 hrs with a standard deviation of ≈4 hrs. It is possible to successfully follow the build-up to eruption in the simulations for all four eruptions (100%), that were observed to originate from low in the corona along the internal PIL by Yardley et al. (2018a).
This indicates that by applying the method of Mackay, Green and van Ballegooijen (2011) to construct a time series of NLFF fields, using the simplest initial and boundary conditions, it is possible to capture the key features of the observable coronal structures in the sample of ARs. To improve on these results the effect on the simulated coronal magnetic field of additional physical effects as well as varying the initial and boundary conditions are examined in the following section.

Consequences of Additional Physical Effects
Although it is possible to simulate the coronal field evolution of an AR using only the LoS magnetic field as the lower boundary condition combined with a potential field as the initial condition, such a simple model does not work in all cases. There were several issues that were encountered in the simulation when using the simplest initial and boundary conditions (an initial potential field condition and closed top boundary). Firstly, the presence of highly twisted field near the side boundaries of the box. Boundary effects can be rectified by re-scaling the magnetograms to occupy a smaller area at the bottom of the computational box during the clean-up procedure (Appendix A). If the magnetograms contain large amounts of small-scale magnetic field that affect the simulated coronal evolution, these can be removed by smoothing the magnetograms with a Gaussian kernel (see Appendix B). This process is applied in addition to the clean-up procedure detailed in Appendix A. If the simulation runs for long time periods, twisted magnetic field can build-up in the computational volume. By adding coronal diffusion, in the form of Ohmic diffusion or hyperdiffusion, this can help prevent the build-up of highly twisted field by decreasing the amount of poloidal flux. However, despite the inclusion of additional coronal diffusion terms, flux ropes are still able to form and reach instability in the simulation and the overall evolution of the simulated coronal field remains significantly unaffected (Mackay and van Ballegooijen, 2006a;Yardley, Mackay and Green, 2018b).
The energy and non-potentiality of the coronal field in the simplest simulation setup only originates from the Poynting flux due to horizontal motions. For the cases where the simple model is insufficient to describe the observations (ARs with a score of 2 or 3) there could be additional physical effects that are acting. For example, the initial configuration of the coronal magnetic field could be non-potential and therefore a LFF field initial condition could be implemented to represent any non-potential effects present before the start of the simulation. When a LFF field initial condition is used the force-free parameter α is assigned a small value with a magnitude of 10 −9 -10 −8 m −1 (see Table 2), to match the weak shear seen in the coronal observations. The range in the force-free parameter is constrained by the size of the computational domain which scales as 1/L, where L varies from one AR to the next. This is due to the nature of the LFF field solution requiring a decaying (non-oscillatory) solution with height. The sign of α is taken from the sense of twist from the magnetic tongues present in the observations (Luoni et al., 2011). The sign and value of α in our simulations is therefore selected in a similar manner to our previous study (Yardley, Mackay and Green, 2018b).
There may also be other sources of energy or helicity injection, which are not captured by the evolution of the normal component of the magnetic field that have to be taken into account, such as the presence of vertical motions or torsional Alfvén waves. Along with these additional injection mechanisms nonideal processes may also have to be considered. These effects are implemented one at a time in the simulation by modifying the induction equation to include the physical effects through three additional terms: The first additional term is Ohmic diffusion, where η represents the resistive coefficient. The second additional term, is hyperdiffusion (Boozer, 1986;Strauss, 1988;Bhattacharjee and Yuan, 1995). This diffusion term is artificial and is introduced to reduce gradients that are present in the force-free parameter α, while total magnetic helicity remains conserved (van Ballegooijen and . The third additional term represents the injection of a horizontal magnetic field or twist component at the photospheric boundary. In this term ∇ z is the vertical component of the gradient operator and ζ is an injection parameter that has the dimensions of a diffusivity. The parameter ζ is only non-zero at the photospheric boundary (z = 0) hence, the injection of the horizontal field only occurs at this location. This term leads to a change in A z half a grid point into the domain, and the subsequent injection of a horizontal magnetic field and magnetic helicity into the corona. By applying this injection in A z leaves the vertical component of the magnetic field unchanged. The sign of the injection parameter ζ determines the sign of the magnetic helicity that is injected via the horizontal field. A positive (negative) value of ζ leads to the injection of negative (positive) magnetic helicity. Once injected, the horizontal field and twist component propagate upwards along the magnetic field lines through the v × B term in the induction equation above (Equation 3). This term is mathematically equivalent to that used in Mackay, DeVore and Antiochos (2014) to model the helicity condensation process of Antiochos (2013). For the present simulations this term does not represent helicity condensation rather it is used to add an additional non-potential contribution that is not captured by a potential field initial condition or the applied horizontal motions on the photospheric surface alone. Additional sources of helicity may originate from the prior evolution of an AR that is not captured from the initial potential field, the presence of vertical motions or the propagation of torsional Alfvén waves from below the photosphere into the corona. The additional injection of horizontal magnetic field at the photosphere, along with the Ohmic and hyperdiffusion terms are included in the simulation through user-defined constants.

Magnetic Field Evolution
We now modify the top boundary, initial condition and include non-ideal terms in the simulations. This is to determine whether it is possible to improve the showing the evolution of AR 12455. The second row (e-h) shows sample field lines from the simulation run with closed top boundary conditions and an initial potential field i.e. the simplest initial and boundary conditions. The third row (i-l) shows the results when Ohmic diffusion, η is added with a value of 25 km 2 s −1 and small-scale field has been removed. The score given to each simulation is given at the bottom right of panels (h) and (l). The positive (negative) photospheric magnetic field is represented by the red (blue) contours.
simulation results, obtained for the ARs in Section 4.1, where only a reasonable or minimal amount of the coronal structure was captured (ARs assigned a scoring criteria of 2 or below).
To improve the results obtained by using the simplest initial and boundary conditions additional physical effects, Gaussian smoothing, and LFF field initial conditions are used. The simulations that were performed for each AR to improve the previous results are described in the comments (final) column of Table 2. If the performance of the simulation improved, the new score is included in brackets in the Score (fourth) column of Table 2. An example can be seen in Figure 5 where AR 12455 improves from a score of 2 to 1. The original simulation captured the evolution of the large-scale coronal loops in the north of the AR relatively well however, failed to reproduce the large-scale loops present in the south. It also failed to capture the bright core of the AR (see Figure 5 (a)). By removing the small-scale magnetic field at the AR periphery using a Gaussian kernel and then introducing Ohmic diffusion the simulation is able to replicate the small and large-scale coronal structure for the entire AR evolution including the sheared structure present at the start of the AR evolution.
The new results are as follows. The enhanced simulations are able to capture the majority of the coronal structure for 12 ARs, a reasonable amount of struc-ture for five ARs and little or no coronal structure for three ARs (see brackets in Table 3). The new results show that one AR moved from scoring category 3 to 2 and two ARs moved from 2 to 1 indicating there was an overall improvement of 20% when a mix of additional physical effects are included. Therefore, the NLFF field simulation is able to capture a reasonable amount of the structure for 85% of the ARs and only failed to capture the structure for 15% of the ARs from the sample. This is a slight improvement on the previous result, where the simplest initial and boundary conditions were used. The improvement in the results is mainly due to the use of a LFF field initial condition. Although, the application of Gaussian smoothing to remove additional small-scale magnetic field near the AR periphery and the addition of Ohmic diffusion also improved the results. When considering the build-up to eruption in the simulations, no improvement is made on the previous results as the simulations again successfully follows the build-up to eruption for 12 out of the 24 observed eruptions.

Discussion
We have used the method of Mackay, Green and van Ballegooijen (2011) to simulate the full coronal evolution of 20 bipolar ARs, from emergence to decay, using a time series of LoS magnetograms as the lower boundary condition. To reproduce the full coronal evolution of the ARs requires a series of magnetograms that extends over the entire lifespan of each AR.
Numerous clean-up processes (see Appendix A) have been applied to the raw magnetograms including time-averaging, removal of isolated features, removal of low flux values, and flux balancing before carrying out the simulations. The application of these procedures produces a series of cleaned magnetograms with a smooth and continuous evolution of the photospheric magnetic field. By using a series of cleaned magnetograms as the lower boundary condition it is easier to simulate the large-scale coronal magnetic field evolution of the ARs as the inclusion of small-scale magnetic elements and random noise could potentially lead to numerical problems in the simulations.
The method has not yet been tested using vector magnetograms as the lower boundary conditions of the simulation. However, an initial qualitative comparison between the vector components at the simulation boundary to the observed vector data of one AR (AR 11561) in our sample shows a relatively good agreement (see Appendix C for more details). In a follow-up study we will expand on this qualitative comparison between the simulated and observed vector magnetic field components.
Initially, the ARs were simulated using the simplest initial and boundary conditions i.e. a potential field initial condition and closed top boundary. We conclude, after a manual comparison with the observations, that the simulations reproduced a reasonable amount of the coronal structure and evolution for 80% of the ARs. This result is improved slightly to 85% by applying Gaussian smoothing to remove additional small-scale magnetic field in the magnetograms, using a LFF field initial condition, and including additional effects such as non-ideal terms in the simulations. For the ARs where the simulation failed to reproduce the main coronal features, particularly during the early stages of the AR evolution, a NLFF field initial condition may be more appropriate. We will implement the use of a NLFF field initial condition in the future by constructing a NLFF field extrapolation using the technique described by Valori, Kliem and Keppens (2005). Whereby a potential field will be extrapolated from a magnetogram, the horizontal field components will be set using a vector magnetogram, and magnetofrictional relaxation will be applied to relax the magnetic field to a force-free equilibrium.
We do not constrain the simulations with coronal observations therefore, the coronal structures reproduced by the simulation are the result of the nonpotential effects produced by the boundary evolution. Therefore, the accuracy of the coronal field models in this study has been judged qualitatively by a manual inspection and visual comparison to the coronal observations. To make the comparison to observations less time-consuming and to remove the subjective nature of this analysis an optimization method could be developed to minimize the deviation between the field lines from the simulation and the intensity observations. An optimization technique will be considered in future studies.
We now discuss when and where the NLFF field simulations were able to reproduce the build-up to eruption. By reproducing the build-up to eruption, we are referring to the ability to identify a flux rope that has formed in the simulation that loses equilibrium or becomes unstable at the same location and at a similar time to the eruption that occurred in the observations. We do not aim to recreate the full dynamics of the eruptions as this requires a MHD simulation.
The simulations were able to replicate the formation and eruption onset of flux rope structures at the internal PIL of an AR where the flux rope was created by flux cancellation and magnetic reconnection occurring at low atmospheric heights. Signatures of eruption onset were found in the simulations for all four low corona eruptions that originated from the internal PIL of ARs 11437, 11561, 11680, and 12382. The simulations were analysed within a ± 12 hr window of the eruption occurring in the observations and the mean unsigned time difference of eruption onset taking place in the simulations compared to the observed eruptions in these four ARs was found to be ≈5 hrs. These simulation results support the van Ballegooijen and Martens (1989) scenario and show that the physical processes can be replicated on a similar timescale to that which the Sun evolves over.
The technique failed to capture the onset of some of the eruptions that originated from low in the corona along an external PIL or at high-altitudes. There are a number of possible reasons for this. Capturing the initiation of eruptions that occurred during the early stages of the simulations proved challenging since these eruptions occur during the flux emergence phase of the ARs. Using an initial potential field condition, combined with the short time over which the coronal field is being evolved, means that insufficient shear and free energy will have built-up in the simulated coronal field.
To combat this issue we can vary the initial or boundary conditions and include additional non-ideal effects in the simulation. For six of the ARs we constructed a LFF field initial condition to see how this affected the results. We chose the magnitude and sign of the force-free parameter to reflect the weak shear seen in the coronal observations. Ideally, vector data can be used to calculate the value of α to use to construct the LFF field initial conditions for the simulation. In the future, we aim to use the observed α value for the LFF field initial condition in our simulations or use a NLFF field initial condition when possible.
The simulation method also fails to capture the eruption onset for ejections that occur in quick succession as it is impossible to separate them from one another in the simulation. To recreate the dynamics of multiple eruptions over short timescales requires the use of full MHD simulations. For example, four eruptions from external PILs were observed to occur in quick succession during the first 12 hrs of the emergence phase of AR 12229. The build-up to these eruptions was not captured by the simulation and this AR accounted for a large number of the missed eruptions. There was also a large imbalance in the magnetic flux during emergence due to the AR emerging at ≈50 • longitude into negative quiet Sun magnetic field. Additionally, eruptions that originated along an external PIL were observationally found to occur due to flux cancellation that takes place between the periphery of the AR and quiet Sun magnetic field during the emergence phase (Yardley et al., 2018a). Eruptions that form at the external PIL are harder to simulate because much of the small-scale field is removed during the "cleaning" process or is not included in the local simulations. At present the simulation method is designed to capture the local and internal evolution of the ARs.
In Yardley et al. (2018a) the origin of each high-altitude event was not studied in detail but it was suggested that they could be the result of the formation of a high-altitude structure during the evolution of the AR or the destabilization of a pre-existing external structure. The build-up to the high-altitude eruptions that were observed in ARs 11446, 11886, and 12336 were not replicated by the simulations. This could be taken into account in future work by using a NLFF field initial condition on a case-by-case basis if a flux rope is present at the start or early stages of the simulation. However, if the high-altitude events are a result of the destabilisation of pre-existing structures then this technique will not be able to capture their formation. Therefore, to be able to capture the onset of eruptions that arise due to the interaction of external magnetic fields or largescale coronal structures, non-local effects need to be taken into account by using global NLFF field models (e.g. Mackay and van Ballegooijen 2006a) to simulate the evolution of the large-scale corona.
Presently, we have focused on simulating the evolution of a set of relatively simple, bipolar ARs that produce faint eruption signatures and a limited number of CMEs. However, this is necessary to test the method before simulating larger, more complex ARs. In the future we will simulate a broader range of ARs, including multipolar regions and large AR complexes that produce multiple CMEs. Given the results of the applied technique simulating larger, multipolar and non-isolated regions should be possible but will require a larger computational domain.

Summary
In this study, the coronal evolution of 20 bipolar ARs was simulated from emergence to decay. The simulations were carried out in order to test whether the evolution of the coronal magnetic field through a series of NLFF states driven by boundary motions could successfully reproduce the observed coronal features of the ARs and the onset of eruption. The coronal magnetic field evolution was simulated by applying the NLFF field method of Mackay, Green and van Ballegooijen (2011) to LoS magnetograms taken by SDO/HMI that were used as the lower boundary conditions. The simulated coronal field evolution for each AR was manually compared to the 171 and 193Å emission structures as seen by SDO/AIA.
The first simulation results were obtained using the simplest initial and boundary conditions i.e. a potential field initial condition and a closed top boundary. By using this approach it was possible to reproduce a reasonable amount of the coronal structure and evolution for 80% of the AR sample. In total, the build-up to eruption was successfully followed in the simulations within a ± 12 hr window of the eruptions occurring in the observations for 12 out of the 24 (50%) of the observed eruptions.
To improve the simulation results we varied the boundary (from closed to open) and initial condition (from potential to LFF) and included additional parameters such as Ohmic diffusion, hyperdiffusion, and an additional injection of horizontal magnetic field and magnetic helicity in the simulations. We also took into account boundary effects by re-scaling the magnetogram at the bottom of the computational box and removed small-scale magnetic features that affect the large-scale evolution of the coronal field by applying Gaussian smoothing to the magnetograms. These steps were in addition to the clean-up procedures and were carried out one at a time. Through considering various combinations of additional terms there was a slight improvement in the results, as one AR moved from scoring category 3 to 2 and two ARs moved from category 2 to 1. Therefore, by varying the boundary and initial conditions and including additional physical effects in the simulation there was an overall improvement of 20%. Overall, the simulations were able to capture a reasonable amount of coronal structure for 85% of the AR sample, only failing to capture the structure of 15% of the regions. Despite varying the boundary and initial conditions and including additional global parameters the simulations are only able to successfully follow the buildup to eruption for 50% of the observed eruptions associated with the AR sample. For the successful cases, the key component in reproducing the coronal evolution and build-up to eruption for the ARs is the use of LoS magnetograms, as the lower boundary conditions to the simulations as changing the side/top boundary conditions, initial condition and including additional physical affects had an insignificant effect on the simulated coronal field evolution.
The unsigned mean time difference between the signatures of eruption onset in the simulations compared to the observed eruptions was ≈5 hrs. The simulations were carried out over a time period of roughly 96-120 hrs therefore, a mean time difference between eruption onset occurring in the observations compared to the simulations of ≈5 hrs is a very favourable result (within 3 applied magnetograms). As current space weather forecasting methods can only provide a warning post-eruption and 1-3 days before the arrival of a CME at Earth with an uncertainty of 12 hrs, our results are well within the present time error. Also, as our approach is computationally efficient we can reproduce the coronal magnetic field evolution of ARs over several days within a few hours of computation time on a desktop machine.
In fact, Pagano, Mackay and Yardley (2019a,b) have demonstrated how eruption metrics based on the NLFF field simulations may be used to distinguish eruptive from non-eruptive ARs. This work has also demonstrated how it is possible to provide near-real time alerts of eruptions using the observed LoS magnetograms, the NLFF field simulations and the projection of the simulations forward in time. The analysis carried out in these studies includes four ARs taken from our AR sample in this paper. The initial results from Pagano, Mackay and Yardley (2019a,b) are promising but additional work is required, including addressing the issues outlined in Section 5, before the method can identify the exact eruption time and be implemented for CME forecasting purposes.
In summary, the full coronal magnetic field evolution of 20 bipolar ARs was simulated using the time-dependent NLFF field method of Mackay, Green and van Ballegooijen (2011). Using this method, it was possible to reproduce the main coronal features present in the observations for 85% of the AR sample. The simulations were also able to successfully follow the build-up to and onset of eruption within a ±12 hr window for 12 out of the 24 eruptions (50%) that were identified in the observations. The mean unsigned time difference between the eruptions occurring in the observations compared to the time of eruption onset in the simulations was found to be ≈5 hrs. It is important to acknowledge that for all four eruptions that took place along the internal PIL of the ARs, the simulations were able to model the timings of eruption onset with a mean unsigned time difference of ≈7 hrs. Therefore, the simulations were able to successfully reproduce the local evolution for the majority of the ARs in the sample.
Acknowledgments The authors would like to thank SDO/HMI and AIA consortia for the data, and also for being able to browse this data through JHelioviewer (http://jhelioviewer.org, Müller et al. 2017). The analysis in this paper has made use of SunPy

Disclosure of Potential Conflicts of Interest
The authors declare that they have no conflicts of interest. Table 1.: The 20 bipolar ARs simulated in this study. The table includes the NOAA number assigned to the AR and the heliographic coordinates of the AR at the time of emergence. The value of peak unsigned flux (half the total absolute positive and negative flux), the start of emergence, peak unsigned flux and end of observation times are also given. The timings of the events that originate from low altitude along the internal PIL, along external PILs, and from high altitude are listed in the final columns. The time and GOES class of four flares and the timings of the two CMEs that are observed in LASCO/C2 that are associated with the ARs are given in the footnotes. The AR properties in this  Table 2.: Results of the NLFF field simulations. The table includes the NOAA AR number, the number of time steps or magnetograms used to simulate the coronal evolution, and whether there was a flux imbalance present between the positive and negative photospheric polarities of the AR. This is followed by the agreement between the simulations and observations, the number of observed eruptions the simulation can follow the build-up to eruption for, and the time difference between the signatures of flux rope eruption onset that occurred in the simulations and the eruptions in the observations. The final column gives additional information such as the location of the AR during emergence if close to 60 • , as well as the surrounding magnetic field environment the region emerges into. It also gives the initial conditions, boundary conditions and additional global parameters that were used to improve the performance of the simulation. If improvements were made, the new score is given in brackets in column four.  No. of ARs 10 (12) 6 (5) 4 (3) Percentage (%) 50 (60) 30 (25)  where C i is the ith cleaned frame and takes values between 1 to n where n is the number of magnetograms in the sequence. F J is the jth raw frame, and τ represents the frame separation where the weighting decreases by 1/e. In this study, the frame separation is set to two meaning that each cleaned frame is a linear combination of the total number of raw frames where the two frames before and after the current frame are weighted the highest. This procedure removes random noise and retains the large-scale features of the ARs. As previously stated, this study focuses on the large-scale evolution of the AR magnetic field and not small-scale elements of the quiet Sun. The next step in the clean-up process includes the removal of small-scale isolated field pixel-by-pixel by evaluating the eight nearest neighbours of each pixel. When fewer than four of the neighbouring pixels have the same sign of magnetic flux then the value of magnetic flux of that pixel is set to zero. Therefore, the pixels at the edge of the magnetogram also have their values set to zero as they have less than four nearest neighbours. In addition, any pixels that have a magnetic flux value below a 25 Mx cm −2 threshold are part of the background magnetic field of the quiet Sun and are also set to zero. At this point the user can choose how to place the magnetograms within the box i.e. the magnetograms can be scaled up/down to fit the computational box or a custom scaling can be applied. In this study, to avoid boundary effects, we rescale the magnetograms to fill 60-70% of the computational box.

NOAA
The last clean-up process is implemented when the top boundary condition in the simulation is set to closed and the magnetograms need to be flux balanced. To flux balance the magnetograms the signed magnetic flux of each frame is calculated. The pixels of non-zero value are summed for each frame and the signed magnetic flux is divided by this total. From every pixel that has a nonzero value the imbalanced magnetic flux per pixel is deducted. As the maximum correction is less than 25 Mx cm −2 no pixels change sign during the balancing Figure 7. The raw and cleaned magnetograms for simulation frame 35 taken on 2013 October 13 when AR 11867 is at its maximum unsigned magnetic flux value. A Gaussian kernel is used to remove small-scale magnetic field surrounding the AR before the clean-up processes described in Section A are applied. For both magnetograms the saturation levels of the photospheric magnetic field are ±100 G. The flux-weighted central coordinates for the positive and negative photospheric magnetic polarities are represented by the red and green asterisks, respectively. of magnetic flux. This is the same threshold that is used to set pixels that form part of the background quiet Sun magnetic field to zero.

B. Gaussian Smoothing
In some cases the small-scale magnetic field at the periphery of the AR is removed before the clean-up procedure is applied (Figure 7). This is achieved by using a method similar to Yardley et al. (2016). Firstly, a Gaussian filter is applied with a standard deviation (width) of 7 pixel units to smooth the data. Secondly, the weighted average of the value of magnetic flux density of the neighbouring pixels must exceed a 40 G cut-off. Then, the largest regions that are identified that make up at least 60% of the selected area are kept whereas, smaller features at large distances from the AR are discarded. This procedure removes small-scale quiet Sun features that are not part of the AR, and does not affect the coronal evolution as we are only interested in the large-scale coronal evolution of the AR.

C. Comparison to Vector Magnetic Field Observations
The method used to simulate the coronal evolution of our AR sample uses a time series of LoS magnetograms as the photospheric boundary condition. This boundary condition injects electric currents into the coronal magnetic field which then evolves through a time series of NLFF fields using the magnetofrictional relaxation process. At no point in the simulation do we constrain the solution using the observed vector magnetic field or with coronal observations. Therefore we allow the boundary evolution to self-consistently produce the horizontal field and subsequently the coronal structures. The reproduced coronal structures are therefore due to non-potential effects produced by the horizontal evolution of the LoS magnetic fields along with any flux emergence or cancellation.
To show that our magnetic field at the photosphere is consistent with the observations and that the simulated coronal structures can be compared with the observed ones we have included a comparison of our simulated vector field at the boundary with the observed vector field. Figure 8 shows the vector magnetic field components from the simulation on the base compared to the observed vector data for AR 11561 where the comparison is carried out midway through its evolution. To produce the observed vector field components we have used the space weather HMI active region patches (SHARPS, Bobra et al. 2014) data that have been projected to the Lambert cylindrical equal-area (CEA) Cartesian coordinate system i.e. the hmi.sharp cea 720s series. The figure shows that there is a relatively good agreement between the simulated horizontal field components and the observed horizontal field, particularly in the strong field regions where the signal-to-noise ratio is high. To determine quantitatively whether there is a close correspondence between the horizontal components derived from the base of the simulation and the observed vector components we will conduct an indepth comparison in a follow-up study. This will include a detailed comparison of the sign and distribution of the three magnetic field components for both the vector field of the simulation and observed vector data. The vertical component of the magnetic field is shown where the positive (negative) photospheric polarities of the AR are represented by the black (white) contours saturated at 500 G. The red arrows represent the magnitude and direction of the horizontal magnetic field components.