© P. Pagano et al., Published by EDP Sciences 2018

1 Introduction

The solar corona is a highly dynamic environment where magnetic and plasma structures are continually evolving. It is a key region of the solar atmosphere that couples the solar photosphere with interplanetary space. This coupling results in the solar wind and many other perturbations of the interplanetary space plasma, referred to as space weather. The origin of these events can be traced back to flux emergence and flows at the solar photosphere. The corona stores vast amounts of free magnetic energy and the release of this energy in the form of coronal mass ejections (CMEs) initiates the largest and most violent perturbations of the near Earth environment, referred to as space weather or extreme space weather (Hapgood, 2011). This topic has drawn the attention of national and international institutions due to the practical consequences that severe space weather events can have on national infrastructures as well as on space operations and missions (Schrijver et al., 2015).

CMEs produce major disruptions to the ambient magnetic and plasma properties of interplanetary space, as their occurrence signifies the violent ejection of both magnetic flux and plasma from the low solar corona. Often magnetic flux ropes in the low corona are considered the main progenitors of CMEs, where an initial period of stability of the flux rope is followed by a fast and sudden ejection outwards into interplanetary space. This scenario is supported by a number of observations (Howard & DeForest, 2014; Chintzoglou et al., 2015; Cheng et al., 2011; Li & Zhang, 2013). Some models explain the flux rope formation and the subsequent ejection with photospheric flows (Mackay & van Ballegooijen, 2006a; Xia et al., 2014), others focus on magnetic flux emergence from underneath the photosphere (Archontis & Hood, 2012), and finally others rely on the onset of MHD instabilities (Török & Kliem, 2005; Zuccarello et al., 2015). It should be noted however, that it is not unusual to observe CMEs bearing no apparent connection with structures in the low corona (D’Huys et al., 2014; Ouyang et al., 2015). A statistical study carried out by Hutton & Morgan (2015) showed that not all CMEs can be directly associated to magnetic flux ropes and Vourlidas et al. (2013) measured that at least 40% of CMEs are certainly associated with magnetic flux ropes, while only a minor percentage are certainly not associated with any magnetic flux ropes. The above discussion illustrates that there is a strong relationship between CMEs and flux ropes.

While the CME initiation mechanism is still under debate, it is widely accepted that, no matter the origin of the CME, they travel outwards due to their own propulsion until about ~4 R_⊙ (Gopalswamy et al., 2000), beyond which they are dragged out by the solar wind. Sachdeva et al. (2015) claim that the aerodynamic drag only becomes the dominant force after 15 R_⊙, although it can play a role at shorter radial distances. It therefore seems reasonable to assume that CMEs are not significantly coupled to the solar wind until after 4 R_⊙. Nevertheless, below 4 R_⊙ and also at larger distances, CMEs can suffer both deflection or acceleration (Gopalswamy et al., 2003). This may be due to the background corona, the interplanetary magnetic field or interaction with other eruptive events.

To mitigate the effects of space weather it is key to predict the arrival time and properties of CMEs at the Earth’s magnetosphere. The development of these predicting tools is not trivial and while current attempts have made significant advances, it is commonly accepted that more precision and more accuracy are required for effective predictions. The most recent development of the drag-based model (DBM) (Vršnak et al., 2013; Žic et al., 2015) is used to predict the arrival time of CMEs assuming that beyond 20 R_⊙ the aerodynamic drag force dominates the dynamics. Moreover, the model ElEvoHI (Rollett et al., 2016) empowers the DBM with an improved geometrical fitting of the CME. Tóth et al. (2005) introduced the Space Weather Modelling Framework (SWMF), a series of highly complex computer simulations that account for significantly different physical regimes and to this end utilises a number of communicating codes (Tóth et al., 2012). The flexibility of the SWMF allows for the coupling of specific models with different codes to develop new modelling techniques, as been done by Jin et al. (2017). Merkin et al. (2016) used two different MHD codes to model the corona and the interplanetary space and a common spherical shell at 20 R_⊙ is used to couple the two codes. Finally, the 3D MHD WSA-ENLIL model describes how the solar wind and interplanetary plasma act as a background for a kinematically inserted CME (Odstrcil et al., 1996; Odstrcil & Pizzo, 1999a, 1999b; Odstrcil, 2003; Odstrcil et al., 2004). In particular, the so called CME Cone Model (Xie et al., 2004; Xue et al., 2005; Michalek, 2006) is widely used in WSA-ENLIL to model the insertion of a CME into the solar wind. Its parameters can be tuned by coronagraph images (Millward et al., 2013), or data collected in situ at Mercury (Baker et al., 2013). Recently it has also been extended to describe Halo CMEs (Na et al., 2017).

For all present and future models the coupling between the outer corona and the lower corona is key. This includes the coupling of the solar wind streams with the origin of CMEs. In particular realistic boundary conditions for space weather need to consider the origin of CMEs in the low corona and how to insert these CMEs into the solar wind. The cone model performs this key role in the WSA-ENLIL model and has improved its predictive capabilities (Dewey et al., 2015), however it injects only a density and velocity perturbation in the solar wind, neglecting the magnetic flux that is an essential part of the CME disturbance to space weather. As an alternative to this, Shiota and Kataoka (2016) introduced the use of a spheromak to model the perturbation from an ejecting magnetic flux rope. It considers the injection of an idealised structure described in terms of density and magnetic field. In Merkin et al. (2016), where two different MHD codes are coupled, the common domain between the two codes is in fact a time dependent boundary condition for space weather. Other techniques to input CMEs into space weather models use coronagraph data directly, such as in the Tappin-Howard (T-H) model (Howard & Tappin, 2009; Tappin & Howard, 2009). In this model ICMEs are reconstructed from visible light images. As an alternative, in Bisi et al. (2013) remote-sensing radio observatios of Interplanetary Scintillation (IPS) are used to probe the inner heliosphere and to identify density irregularities in the solar wind. When enhanced by the UCSD time-dependent tomography technique (Jackson et al., 2010). Harrison et al. (2017) provides an extensive review on how heliospheric imaging can be used to improve our space weather forecasts.

Many recent studies have measured our predictive capability. Zhao & Dryer (2014) estimated the uncertainty in our capability to predict the arrival time of CME at 12 hours. Also Falkenberg et al. (2011) found a similar result for WSA-ENLIL. They explain that one of the reasons that prevents better estimates is our limited capability in reproducing solar transients and how the CMEs are input in the background solar wind. This seems a general trend, as no major difference has been found between the various models, where, for instance, the DBM and ENLIL predicting capabilities differ by less than 10%, except for the case of strong solar activity when ENLIL performs significantly better (Vršnak et al., 2014). More recently Tucker-Hood et al. (2015) carried out an extensive survey that proved how much room for improvement there is in the field. Out of the 60 predictions made, 36 were false alarms and when the prediction was successful the arrival time was still estimated with an error of ~16 h. Mays et al. (2015) came to similar conclusions, but found that results were slightly more positive for ENLIL. These results show that significant modelling improvements are needed to reduce the arrival time error and the property of the geomagnetic perturbations. To improve these estimations more accurate boundary conditions of the injection of CMEs into interplanetary space are required.

In the present paper we focus on a specific kind of possible improvement: to provide realistic and accurate boundary conditions for the next generation of space weather forecasting models. In particular we need realistic initial and time dependent boundary conditions that reflect the complexity of the transition from the solar corona to interplanetary space, in terms of the injection of both plasma (density and velocity) and magnetic field in the solar wind. To do this we present a novel approach where simulations of magnetic flux rope ejections that are derived directly from surface magnetograms may be used to aid future space weather predictions. While most magnetic structures in the solar corona are close to equilibrium, it is essential to identify those that lose equilibrium and then erupt. To this end Yeates et al. (2007) have opened the way for the use of synoptic magnetograms to model the quasi-static, non-potential evolution of the global corona. This modelling technique has proved to be accurate in predicting the formation and helicity of solar filaments (Yeates et al., 2008), the variation of the Sun’s open magnetic flux (Yeates et al., 2010b) and to a lesser extent the location of CMEs (Yeates et al., 2010a). As already described in Mackay & van Ballegooijen (2006a) and Pagano et al. (2013b), the formation and ejection of magnetic flux ropes occur over very different time-scales and in different dynamic regimes. The formation of magnetic flux ropes occurs slowly (days or weeks), through a series of quasi-equilibrium states (formation times much longer than the Alfvén time) and in a magnetically-dominated regime (β ≪ 1). In contrast, magnetic flux rope ejections takes place over a few hours, where no equilibrium exists and where compression and heating create extended regions of high-β plasma. To simulate the global corona during both the formation and eruption of magnetic flux ropes, we expand upon the work already carried out in (Pagano et al., 2013a, b, 2014). In these studies we have shown that coupling a quasi-static non-potential model with MHD simulations is a viable way to describe the full life span of a flux rope: from formation to ejection. One limitation of the previous studies, was that they only considered a small wedge-shaped portion of the Sun and simple idealised magnetic field configurations. We now apply the same approach to realistic observationally derived magnetic fields that cover the full sphere of the Sun as simulated in the global non-potential model of Yeates et al. (2010a). The model of Yeates et al. (2010a) uses bipoles deduced from synoptic magnetograms to realistically model both the time evolution of photospheric fields and the quasi-static response of coronal fields to global motions over long periods of time. One aspect of this is the formation of magnetic flux ropes, which subsequently lose stability. At this loss of stability we switch the modelling approach where MPI-ARMVAC is used to follow the MHD evolution. Throughout the MHD simulation we derive the boundary conditions for space weather forecasting tools, i.e. the density, velocity and magnetic field distributions that the flux rope ejections inject in the solar wind and interplanetary space.

In the present paper we present the basic physics of a model that aims to produce boundary conditions for space weather prediction models. At the present time we focus on the properties of the model, where in future studies we will consider its application and viability in the context of an operational model along with its coupling to interplanetary models. The significance of this approach relies on its accuracy and efficiency. The global non-potential model has been extensively tested and has been shown to be accurate in describing the evolution of the photospheric magnetic field and subsequent formation of non-potential magnetic fields and flux ropes in the corona. General MHD simulations (capable of handling multi-β domains) are the most accurate tool to model flux rope ejections and then the evolution of the flux rope through interplanetary space. While we put forward a two-stage approach for producing the boundary condition it is of course possible to use MHD simulations to simulate both the magnetic flux rope formation and ejections. However, present computational power is insufficient to model the slow formation of flux ropes over days to weeks. A MHD simulation would consume unreasonable computational resources and additionally the high number of time steps required to cover the physical time span would lead to the accumulation of significant round-off errors. Taking into account these considerations, we propose a three stage numerical model. For the present paper we consider the first two steps of formation and eruption. The third and final stage which is the connection of the ejection stage to an interplanetary evolution model will be carried out in the future. This will include using the output from the MHD simulation as a time-dependent boundary condition for driving a new generation of space weather forecasting tools.

The structure of the paper is as follows: in Section 2 we present the construction of the model, in Section 3 we describe in detail the MHD simulation, in Section 4 we explain how this study can be used in a space weather forecasting context and we finally discuss and give conclusions in Section 5.

2 Model

2.1 Overview of coupled modelling technique

In order to model the ejection of magnetic flux ropes in the global corona, we employ a dual modelling technique of MHD simulations coupled with a quasi-static non-potential global model (Yeates et al., 2010a). This approach is an extension of the technique successfully pioneered in Pagano et al. (2013b) and then further developed in Pagano et al. (2013a) and Pagano et al. (2014). In these studies, a magnetic configuration obtained from the global non-potential model case study of Mackay & van Ballegooijen (2006a) was used as an initial condition in the MHD simulation. The combined technique allows us to follow the slow build-up of stress and electric currents over observed solar timescales (days -months), along with being able to follow the dynamic eruption timescale (min-hours).

2.1.1 Quasi-static model

In this paper, we use a magnetic configuration obtained from a simulation run of the global non-potential model, where the model simulated the entire time span of Cycle 23 (Yeates & Mackay, 2012). A magnetic configuration is chosen near the end of the run, where the structure and connectivity of the corona was produced by the combined effects of differential rotation, meridional flow, surface magnetic diffusion and magnetic flux emergence. The properties of the flux emergence events were deduced from NSO/KP and NSO/SOLIS magnetograms and were included to maintain the accuracy of the photospheric field compared to that found on the Sun. This application of the Global Model has proved to be very successful in reproducing the chirality of solar filaments (Yeates & Mackay, 2012) throughout the solar cycle.

Additionally, we also took care by selecting an initial condition for the MHD simulation where there were two relatively isolated magnetic flux ropes that had formed and both had accumulated enough stress that they were going to erupt. In the global model we have developed an automated technique for the identification of erupting flux ropes. The technique is described in Yeates & Mackay (2009a). The global magnetic configuration under investigation can be seen from two viewpoints, 180° in longitude apart, in Figure 1. The magnetic field lines exhibit a number of features characteristic of those on the Sun, both at the photosphere and in the magnetic field connectivity in the low corona. At several locations opposite polarity patches lie close to one-another, and these are the preferred locations of flux cancellation where flux ropes are found. However, at only two of these locations (FR1 in the northern hemisphere and FR2 in the southern hemisphere) is the magnetic field sheared enough along the polarity inversion line to show a twisted magnetic flux rope. At FR1 a flux rope sits above a PIL in a bipole that extends over a region of . The PIL is mainly North-South directed. At the location FR2 there is a smaller flux rope and bipole, where now the bipole extends only for . The polarities of the bipole are separated by a PIL which lies directed in a South-West to North-East direction. For both of the active regions the flux rope lie above the PIL therefore have a slightly different orientation. The rest of the coronal magnetic field is mostly close to potential and the other bipoles that are present in the domain show no significant shear above their polarity inversion lines.

Fig. 1 Representation of the magnetic configuration imported from the quasi-static non-potential evolution model. The two panels show maps of the radial component of the magnetic field on the surface from two opposite points of view. Blue lines are some representative magnetic field lines to highlight existing structures.

2.1.2 Coupling to MHD simulation

To use the magnetic configuration illustrated in Figure 1 as the initial condition for the MHD simulation, we import into the MHD simulation all three components of the magnetic field from the global non-potential model. A full description of how this is carried out, where both the stability (or instability) and magnetic connectivity of the configuration is preserved, is given in Pagano et al. (2013b) and Pagano et al. (2013a) where a number of test cases are described.

For the present paper we have slightly modified the way in which the 3D interpolation is performed. In spherical coordinates, let A(r, θ, ϕ) be the value of a function that we want to have interpolated to the position (r, θ, ϕ), where r is the radial distance from the centre of the Sun, θ is the polar angle and ϕ is the azimuthal angle. We know that it lies in the cell defined by the indexes [i : i + 1, j : j + 1, k : k + 1] where the original function a is defined, we compute A(r, θ, ϕ) as (1)

where V [i, j, k] is the volume defined by the point (r, θ, ϕ) and the cell corner opposite to the position [i, j, k] and V is the sum of all the volumes. This approach guarantees the continuity of the interpolated solution and its smoothness, independent of the spatial resolution of the grid where A(r, θ, ϕ) is defined.

2.2 MHD Simulation

2.2.1 Plasma distribution

As the global non-potential model only provides a magnetic configuration, to produce a complete set of MHD variables in the initial condition we need to determine distributions of plasma density, velocity and temperature. In specifying these we aim at a realistic and general representation of the solar corona, where the distribution of plasma takes into account the heterogeneity of the structures of the solar corona. This heterogeneity can be simplified by representing the solar corona as a combination of dense active regions and less dense quiet Sun regions. In addition, there are dominantly horizontal magnetic fields (such as prominences or filaments) which are generally two orders of magnitude cooler and denser than typical coronal values.

We define the following proxy function to link the plasma temperature and density to the magnetic field: (2)

where (3) The function ω is positive definite and peaks where the magnetic field exhibits a complex twisted field, e.g. near the axis of a magnetic flux rope. As the value of ω is also proportional to the magnetic field intensity, it is higher near the solar surface (where the magnetic field is more intense) and lower at further radial distances from the solar surface.

In order to effectively use ω to model the solar atmosphere we define the functions: (4) (5) (6) where Ω_θ and Ω_B are functions bound between 0 and 1 while Ω is defined to pick the higher between Ω_θ and Ω_B. Ω_θ is a function that only depends on the coordinate θ that defines two regions near the poles where we quench any dynamics close to the boundaries of the simulation. Using Ω the temperature is defined by: (7)

Next, the thermal pressure is independently specified by the solution for hydro-static equilibrium with a uniform temperature set equal to T_corona, (8)

where ρ_LB is the density at r = R_⊙ when |B| = 0, μ = 1.31 is the average particle mass in the solar corona, m_p is the proton mass and k_B is Boltzmann constant. Finally, the density is simply given by the equation of state applied to Eqs. 7 and Eq. 8: (9)

This produces an inhomogeneous solar corona of cool dense flux ropes and hotter emptier corona arcades.

2.2.2 MHD simulation

Using the approach described in Section 2.2.1, we construct the initial condition for the MHD simulation, where Table 1 shows the value used in our model for all parameters.

With this initial condition we use the MPI-AMRVAC software (Porth et al., 2014), to solve the MHD equations where external gravity is included as a source term, (10) (11) (12) (13) where t is time, ρ is density, v velocity, p thermal pressure and B the magnetic field. The total energy density e is given by (14)where γ = 5∕3 denotes the ratio of specific heats. The expression for solar gravitational acceleration is given by (15)where G is the gravitational constant, M_⊙ denotes the mass of the Sun, and is the unit vector.

The computational domain is composed of 256 × 256 × 512 cells, distributed on a uniform spherical grid. With this resolution the simulation domain extends over 3 R_⊙ in the radial direction starting from r = R_⊙. The co-latitude, θ, spans from θ = 0.75° to θ = 179.25° and the longitude, ϕ, spans 360°. The boundary conditions are treated with a system of ghost cells and match those used in Yeates et al. (2010a). Open boundary conditions are imposed at the outer boundary, reflective boundary conditions are set at the θ boundaries and the ϕ boundaries are periodic. The reflective θ boundary condition does not allow any plasma or magnetic flux to pass through. The stability at this boundary is reinforced through the use of the function Ω_θ , which results in plasma near the θ boundary having a density significantly higher than that prescribed by gravitational stratification. This enhancement quenches any upward motion at this boundary. At the lower boundary which represents the photosphere, we impose a fixed boundary condition taken from the first two θ-ϕ planes of cells derived from the global non-potential model.

The interpolation technique used to import the magnetic field configuration from the spatial grid of the global non-potential model, to the alternative grid that we use in the MHD simulation only applies to r < 2.5 R_⊙. In the region beyond r = 2.5 R_⊙ we assume a purely radial field where magnetic flux is conserved: (16)

The key properties of the initial condition in the MHD simulation can be seen in Figure 2 which shows in Mollweide projection distributions of (a) the radial component of the Lorentz force at the lower boundary, (b) the function ω derived from the magnetic configuration and finally the resulting distributions of (c) plasma density and (d) temperature. The Mollweide projection has the advantage of showing in a single field of view the whole solar disk and therefore clearly illustrates the connectivity between different regions. The Lorentz force exhibits several areas of strong outward radial force shown as the white regions. In particular two structures show a significantly larger positive Lorentz force compared to surrounding areas. One lies in the northern hemisphere to the left of central meridian and the other is in the southern hemisphere to the right of central meridian. These are at locations where large flux ropes have formed and can no longer be fully contained by the overlying fields. Close to where the flux ropes lie, with their positive radial component of Lorentz force, smaller regions where the Lorentz force is negative can be seen.

The map of ω (Fig. 2b) follows a similar distribution, where several small structures are visible along the active latitudes. Three zones exist where ω is significantly larger with peak values nearly ~100 times greater than ambient background values. Flux ropes show higher values of ω, as close to flux ropes all the three terms ω_r , ω_θ , and ω_ϕ are important at different locations. Namely ω_r peaks at the center of the magnetic flux rope axis, while ω_θ and ω_ϕ peak near the footpoints. Finally, the maps of density (Fig. 2c) and temperature (Fig. 2d) show that the locations with a more complex magnetic configuration are denser and cooler. With this the magnetic flux ropes in the simulation appear as regions with density and temperature one or two orders of magnitude denser and colder than the surrounding ambient values. Thus the initial conditions produce an inhomogeneous corona that exhibits many features found on the Sun.

The initial plasma β in the simulation ranges between β ~ 10⁻³ at the flux ropes to β ~ 1 in confined regions where the magnetic field is weak. While there is a wide range, normally the value of β lies between 0.1 and 0.01 throughout the vast majority of the volume. Due to the low β values at the flux ropes, the strongest unbalanced force in the initial condition is the radially directed Lorentz force at these locations. In addition to the unbalanced Lorentz force the radial profile of density and pressure also does not prescribe a balance between the thermal pressure gradient and gravity. While this pressure force imbalance exists, the consequences have been discussed in detail in Pagano et al. (2013a), where it is shown that any evolution due to these imbalances occurs over timescales much longer than that of the eruption dynamics triggered by the Lorentz force. As such we can neglect these effects.

Fig. 2 Mollweide projections (central meridian at longitude ϕ = 180∘ ) on the lower boundary of the MHD simulation of (a) Lorentz Force, (b) ω, (c) Log10(ρ[g∕cm3]) and (d) Log10(T[K]).

3 Simulation

The MHD simulation produces two individual eruptions that originate from the two twisted magnetic flux ropes shown in Figures 1 and 2. In each case, the radial Lorentz force excess due to a local non-equilibrium triggers the dynamic evolution of the eruption that occurs in the simulation. Figure 3 shows images of the column density where the central meridian longitude is chosen to be ϕ = 90° (left-hand side column) and ϕ = 270° (right-hand side column) at (a) t = 0 min, (b) t = 33.6 min and (c) t = 60.3 min where magnetic field lines (green) have been drawn from the photospheric boundary.

The two flux ropes are located such that one is in the northern hemisphere, while the other is in the southern hemisphere, and they lie approximately 180° in longitude (ϕ) apart. When the ejections occur they displace dense plasma to larger radii, which results in an increased column density around the locations of each flux rope. This increased density can be seen as the two circular structures expanding at either limb, one towards North and ϕ = 0°, the other towards South and ϕ = 180°. During the ejections the global magnetic field undergoes a rapid evolution. Even though there is a rapid evolution, it is possible to follow the evolution of the twisted magnetic field lines of the flux ropes. After t = 60.3 min the two circular expanding structures of density have expanded sufficiently that they cross the outer boundary of the domain.

To better understand the evolution of the two flux rope ejections we plot radial cuts of density, temperature, radial velocity and ω from the photospheric boundary to the outer boundary (Figs. 4 and 5). These plots are taken along a radial line intersecting with the initial position of the centers of the two magnetic flux ropes.

Figure 4a shows the density at t = 0 min (black line), t = 10.4 min (blue line) and t = 27.8 min (red line) for both Flux Rope 1 and 2. In each plot the black line shows the initial density profile with the dense flux rope at low heights along with the stratified atmosphere. At later times the presence of the ejection is identified by an excess in density compared to that of the initial stratified atmosphere. At t = 10.4 min the density excess of the ejected plasma has reached, ~1.8 R_⊙ for FR1 and ~2 R_⊙ for FR2 showing that they are traveling at different speeds. Similar estimations can be made with the fronts of the two ejections and at different times. In all cases, we consistently find that FR2 travels about 25% faster than FR1.

At the lead front of the density excess there is a very sharp temperature increase to T ~10⁷ K (Fig. 4b). Beneath the front the temperature profile also shows significant structure. At t = 10.4 min the temperature profile in the atmosphere exhibits a dip at the location of the axis for both FR1 (1.2 R_⊙) and FR2 (1.4 R_⊙). However, it should be noted that the temperature at the center of the flux rope is comparable to that of the background external corona at 2MK. While this value is high, it is still cooler than that of the surrounding ejection.

At t = 10.4 min the radial velocity profiles in Figure 5a for both FR1 and FR2 exhibit a very similar behavior, where there is an initial sharp rise at the photosphere, followed by a slight dip before becoming roughly level between 1.5–2 R_⊙. In both cases at the lead edge there is a sharp fall to zero. Combining the information in Figures 4 and 5 it can be seen that the ejection of both flux ropes can be characterized by a front (density excess, sharp increase in temperature and radial velocity) and a core (density excess, temperature dip, slower than the front). While this is characteristic of the early stages of the ejection similar structures are visible in the plots at t = 27.8 min. However, there are some differences from the earlier times which include: i) the volume of density excess increases and thus the density at the flux rope decreases, ii) the temperature at the front increases, while the temperature dips increase in number and become colder and iii) the radial velocity at the front slightly increases and with this the speed ofthe flux rope ejection. In addition to the features above, we also find that for both FR1 and FR2 ω shows an increase in value with respect to that of the initial condition. This occurs both in the region beyond the front but more prominently near the location of the flux rope center.

The occurrence of the two ejections has a significant effect on the energy budget of the solar corona. Figure 6 shows the total energy as a function of time in the MHD simulation (black line). The total energy is approximately conserved and after t = 68.44 min the total energy has increased by less than 9%. The increase can be accounted for as there is a net flux of energy into the simulation after taking into account the energy flux across the lower and outer boundaries. At the start, the magnetic energy accounts for about 90% of the total energy and the thermal energy for the remaining 10%. By the end of the simulation the kinetic energy is 6% of the total energy where its increase is due to the conversion of magnetic energy which falls to 80% of the total energy, while the thermal energy slightly increases.

Fig. 3 Maps of the column density seen from either side of the Sun in the MHD simulation along with some superimposed magnetic field lines traced from the lower boundary (green) at t = 0 min (a), t = 33.6 min (b), and t = 60.3 min (c). Central meridian for left-hand side panels is at ϕ = 90∘ and at ϕ = 270° for right-hand sided panels and the FOVs are chosen to be same as in Figure 1.

Fig. 4 (a) Cuts of Log10(ρ[g∕cm3]) along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line). (b) Cuts of Log10(T[K]) along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line).

Fig. 5 (a) Cuts of radial velocity along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line). (b) Cuts of Log10(ω) along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line).

Fig. 6 Energy in the MHD simulation integrated over the whole spatial domain as a function of time. Black line shows total energy, blue magnetic energy, red kinetic energy, and green thermal energy.

4 Toward a space weather application

One of the goals of the present work is to provide accurate boundary conditions of the outer solar corona that can be used in future space weather forecasting tools, such as those for the solar wind and the evolution of Interplanetary Coronal Mass Ejections (ICMEs). To achieve a realistic model of space weather conditions at 1 AU, it is key that accurate initial conditions for the injection of plasma and magnetic flux into the solar wind are specified close to the Sun. These boundary conditions require accurate time-dependent modelling of low coronal magnetic fields during the build up to and occurrence of an eruption.

Our technique of coupling the global non-potential model with the MHD numerical solver MPI-AMRVAC is suitable for space weather forecasting purposes because it accurately models observed magnetic field configurations on the Sun and is computationally efficient. It allows us to model the fast dynamic ejection of magnetic flux ropes, maintaining the accuracy and generality of a full MHD model, while allowing us to model the slow quasi-static formation of flux ropes and the global corona in a numerically cheaper way. Initial estimations suggest that the global non-potential model is between 10⁴ and 10⁵ times faster than the MPI-AMRVAC in advancing the simulation in physical time. This therefore allows for the near real time simulation of observed magnetic fields on the Sun. We now discuss how the dual approach described above may provide useful boundary conditions for space weather models.

4.1 Injection of the CME into interplanetary space

In this section, we present some exploratory results on how our combined model injects a CME into interplanetary space. As we have set the outer boundary condition at 4 R_⊙ we focus on the conditions of the plasma and magnetic field as a consequence of the eruption before it starts to interact with the solar wind.

Figure 7 shows Mollweide projections of (a) density and (b) the radial velocity at the outer boundary of the simulation at a time near when both flux ropes cross this boundary (t = 60.3 min).

The two flux rope ejections appear as density excesses covering an area of some tens of squared degrees where FR1 involves a wider area than FR2. One predominately lies in the northern hemisphere and the other in the southern hemisphere, but both straddle the equator. The density of the flux ropes is around 100 times larger than that of the background corona, which has a density of around ~10⁻¹⁹ g∕cm³. Due to the two orders of magnitude difference it appears that at the outer boundary the flux rope has a near homogeneous density. However, on closer inspection of the data internal variations inthe density are also visible, where higher density structures lie in the centre of the ejection for FR1, and at the southern boundary for the ejection of FR2. A similar pattern can also be seen in the radial velocity maps, where both ejections exhibit an outward velocity in the order of hundreds of km∕s. Outside of the region corresponding to the erupting flux ropes the corona has a significantly smaller outward velocity. It should be noted that in this simulation the fastest speeds are not found at the flux rope locations, but rather along the inversion line at the outer boundary, where it is possible that the eruptions trigger some magnetic reconnection processes that accelerate the plasma. At this location where the radial component of the magnetic field changes sign, the radial outward velocity reaches ~1000 km∕s. This outflow is related to a streamer-like phenomenon, where plasma accelerates in conjunction with reconnection at a current sheet. This is partly due to the magnetic field from the global simulation relaxing at this location in the high corona when the MHD simulation starts. This is an artefact of the coupling of the two codes that we will tackle in the future. Whilst this occurs, the low density associated with this feature means that it does not contribute significantly to the momentum flux and is essentially negligible. The internal structure of the radial velocity patterns at both flux ropes differs. For FR1 there are higher radial velocities at the borders of the region involved in the ejection and the region of high positive radial velocity extends to near the polarity inversion line. In contrast the radial velocity pattern of FR2 shows more modest internal structuring and remains isolated from the polarity inversion line.

Figure 8 shows cuts along the ϕ direction at the outer boundary at three different co-latitudes for (a) density, (b) radial velocity and (c) temperature. In each of the plots different co-latitudes are considered corresponding to the center of FR1 (θ = 71.5°, top), the equator (θ = 90°, middle) and finally the center of FR2 (θ = 110°, bottom). In each plot values are shown at three different times: t = 0 (black line), t = 33.6 min (blue line) and t = 60.3 min (red line). The plots allow us to show more quantitatively the inhomogeneous nature at the outer boundary due to the two ejections. The latter two times correspond approximately to the times when the fronts of the ejections reach the outer boundary and when the centres of each flux rope cross it. It is only approximate as the two fronts and the two centres cross the boundary at slightly different times. The plots of (a) density and (b) radial velocity both show that the perturbation of the quiet corona increases in amplitude as the flux ropes reach and cross the outer boundary. When the front of each ejection first reaches the outer boundary (t = 33.6 min) the density increases by less than an order of magnitude and is localised in space. However once the flux rope crosses the surface (t = 60.3 min) the density increases by two orders of magnitude. In contrast, the radial velocity exhibits a number of different features. When the fronts from both ejections reach the outer boundary a peak of ~1000 km∕s occurs, after which the radial velocity settles to a more plateau structure of around 500 km∕s. Variations in both the density and velocity occur over a wide angular degree. Even though both of the initial ejections occur far from the equator, their spatial extension at the outer boundary is wide enough to cross into the opposite hemisphere. From the cuts taken at the latitudes of either FR1 or FR2 there are clear indications of the occurrence of a CME in the opposite hemisphere. In contrast the equatorial cut shows a similar contribution from both FR1 (ϕ ~ 50° in Fig. 8) and FR2 (ϕ ~ 250°) in Fig. 8), as well as from the polarity inversion line plasma flows (ϕ ~ 140°). Finally the temperature of the plasma crossing the outer boundary qualitatively changes, as first the front crosses the equator and then the flux rope core. As the fronts from both ejections cross the boundary, a temperature increase is found. However as the flux rope cores pass through the boundary this changes to a temperature decrease. This behavior seems reasonable from a qualitative point of view, as we expect to find higher temperatures at the front due to the compression and heating of the plasma and lower temperatures at the flux ropes core due to the cooler denser material that is present within it. While we see these temperature features we should also acknowledge that the present model is not designed to model accurately the temperature evolution as it lacks crucial terms such as thermal conduction and radiative losses in Eq. (13). Even though these terms are not included, they are not essential for the modelling of the thermal structure of the front and core of the flux rope ejections.

Figure 9 shows Mollweide projections of the three components of the magnetic field for (a)–(c) the initial condition of the MHD simulation at the level of the photosphere and (d)–(f) at the outer boundary, and (g)–(i) at the outer boundary at the time when both flux ropes cross the outer surface (t =60.3 min). The first row shows B_r , the middle row B_θ and the bottom row B_ϕ. Upon comparing the individual field components at the different heights and at the different times, it is not surprising that the complexity and structure of the magnetic field at the solar surface (Fig. 9a–c) is not carried out to the outer corona (Fig. 9d–f). It is however essential to notice that the configuration of the magnetic field at the outer boundary is the result of the interaction between the outward travelling structures and the background solar corona. This interaction simplifies the configuration of the magnetic field and alters the surviving structures in a way that is not predictable a-priori. In particular while many magnetic features are visible in Figure 9a only two distinct ones appear at the outer boundary during the eruption (Fig. 9g). At the outer boundary the flux rope PILs are of different shape and size, even though the two flux ropes are initially of the same size. Finally, both ejections show a large negative B_θ patch at the outer boundary, but opposite and different B_ϕ patches showing that they undergo a different rotation.

By comparing Figure 9g–i with Figure 9d–f, it is clear that the magnetic field carried by the flux ropes significantly perturbs the outer boundary. The radial component of the magnetic field exhibits at t = 0 s a clearly defined distribution, where there are two regions with opposite polarity that are separated by a polarity inversion line. Within these locations when the flux rope ejection occurs and reaches the outer boundary it alters the magnetic field distribution. The distribution becomes more in-homogeneous and structured. Within each region that was originally of a single polarity, opposite polarity patches now appear indicating the flux rope ejections. In addition to the changes in the radial field, as the initial magnetic field is prescribed to be only radial at the outer boundary, the θ and ϕ components at the outer boundary show a major change once the erupting flux rope reaches it. Both erupting flux ropes carry a negative θ component creating two areas (green) clearly visible in a more uniform and less intense background. Finally, a clear signature of the flux ropes can also be seen in the ϕ component of the magnetic field, where once again two very visible regions are present. An interesting feature is that around FR1 a distinct ring structure of positive B_ϕ is present, within which there is contained a weak B_ϕ. In contrast around FR2, there is a more uniform negative B_ϕ. The origin of this variation is still under investigation, however it illustrates the importance of modelling realistic photospheric magnetic field configurations both before and during the eruption as variations at the outer boundary must ultimately be due to variations in the low coronal field. Thus within the present simulation due to the low coronal properties of the field and the regions in which the flux ropes form, we find significant differences in the field distribution at the outer boundary. If the low coronal field is not accurately modelled then such a variation in ICME’s will not be taken into account.

In Figure 10 we focus on the time evolution of B_θ at the latitudes of FR1 (top), the equator (middle) and FR2 (bottom). Figure 10a shows the quantity at the level of the photosphere at time t = 0.0 just before the eruption starts. Figure 10b shows the results for the outer boundary at the times of t = 0.0 (black line), t = 33.6 (blue line) and t = 60.3 (red line). Again comparisons between Figure 10a,b show no evident correlation between the B_θ values at the lower boundary before the onset of ejections to that of the resulting configuration at the outer boundary as the eruption passes through it.

In Figure 10b the θ component of the magnetic field is particularly important as it gives the out-of-ecliptic plane component at the equator. For space weather consequences the intensity and orientation of this component when the ICME encounters the Earth magnetosphere determines the geo-effectiveness of the disturbance. At the present time it is not clear if the out-of-ecliptic component of the magnetic field at 4 R_⊙ is preserved during the transit of the CME to Earth. However in order to produce accurate predictions at Earth of its intensity and orientation accurate initial conditions are required when it is ejected into interplanetary space. In Figure 10b we find that the ejection front generates a modest enhancement in the B_θ component at all three latitudes. However the overall effect of each front is quite wide. At the equator when the two fronts reach the outer boundary we find a preference for negative B_θ at the longitudes corresponding to the flux ropes. A more significant B_θ enhancement occurs when the flux ropes themselves reach the outer boundary. This is most intense at the center of the flux rope locations, but still visible at other longitudes. At the equator, we find a strong negative B_θ value at the longitudes of the two flux ropes, with a positive peak between them.

Finally in Figure 11 we show a Mollweide projection of ω at the outer boundary at t = 60.3 min when the flux rope crosses the surface. This is to determine whether or not the function ω is appropriate to identify the flux ropes at higher radial distances. Analysis indicates that the locations involved with the flux rope ejections show up as visible features in the ω distribution. These generally occur at high ω values even if it is not possible to define a simple threshold that isolates the flux rope. It is also clear from Figure 11 that thepolarity inversion line is also a very clear feature due to the relatively lower intensity and less twisted configuration.

Fig. 7 Mollweide projections (central meridian at longitude ϕ = 180°) at the outer boundary of the MHD simulation at t = 60.3 min of (a) ρ[g∕cm3] and (b) radial velocity of the plasma [km∕s].

Fig. 8 Cuts of (a) Log10(ρ[g∕cm3]), (b) radial velocity [km/s], and (c) temperature [K] along ϕ at t = 0 min (black line), t = 33.6 min (blue line) and t = 60.3 min (red line) atthree different co-latitudes: the center of FR1 (θ = 71.4°), the equator (θ = 90°), and the center of FR2 (θ = 110°).

Fig. 9 Mollweide projections (central meridian at longitude ϕ = 180°) at the lower boundary of the MHD simulation at t = 0 min (left column) and at the outer boundary at t = 0 min (central column) and t = 60.3 min (right column) for the three magnetic field components Br, Bθ, Bϕ.

Fig. 10 Cuts of Bθ along ϕ at three different latitudes: the center of FR1 (θ = 71.4°), the equator (θ = 90°), and the center of FR2 (θ = 110°). a) At t = 0 and lower boundary. b) At t = 0 min (black line), t = 33.6 min (blue line) and t = 60.3 min (red line) and outer boundary

Fig. 11 Mollweide projection of ω at the outer boundary of the MHD simulation at t = 60.3 min (central meridian at longitude ϕ = 180°) .

5 Discussion, outlook, and conclusions

The work presented here is the first step of a longer-term effort to improve our predictive capability for the occurrence of events on the Sun relevant to space weather. As part of this process we will also improve our understanding of the evolution of non-potential magnetic structures in the solar corona and how they lead to CMEs.

5.1 Present results

In this paper we have carried out global simulations of the Sun where the entire extent of the corona is considered from the photosphere at R_⊙ out to 4 R_⊙. Two models have been coupled to follow the full life cycle of flux ropes from formation to eruption. Initially, a quasi-static non-potential model is used to consider the slow evolution of the coronal magnetic field and formation of flux ropes up until the point of eruption. To follow the dynamics of the eruption an MHD simulation is then carried out. As the quasi-static non-potential simulation does not provide the MHD simulation with a density profile, we put forward a new technique to prescribe the density using the twisted nature of the magnetic field. This allows us to produce an inhomogeneous corona, along with flux ropes that have a density 1-2 orders of magnitude greater than that of the surrounding corona. The simulation presented in this paper shows the ejection of two separate magnetic flux ropes that have formed in the global solar corona at different locations. These flux ropes have been formed self-consistently in magnetic field configurations determined directly from observations. In the simulations the onset of eruption occurs once the tension of the overlying arcades is insufficient to hold down the underlying flux rope. This occurs depending on the size of the flux rope and the strength/topology of the overlying arcades and represents the occurrence of a non-equilibrium. When a non-equilibrium occurs and the flux rope starts to rise, reconnection then occurs under it which results in a radially outward Lorentz force which ejects the flux rope. This can be seen as a two step process of initial slow rise, followed by a faster ejection (see Mackay & van Ballegooijen, 2006a, 2006b for a description).

In the present model the flux ropes are ejected out of the computational domain as the Lorentz forces generated in the corona due to the surface magnetic reconnection and shear are too large to be contained by the overlying arcades (Yeates & Mackay, 2009b). During the ejections, both flux ropes show velocity values along with magnetic and density structures compatible with standard CME parameters at 4 R_⊙.

When the ejections reach the outer boundary at 4 R_⊙ we discuss how the CME passes over this boundary and is injected into the solar wind. We find that the flux rope ejection into the solar wind occurs over an area of tenths of square degree, where the ejected plasma is denser compared to that of the quiet corona. In addition the ejected plasma and magnetic field show a highly inhomogeneous structure. We highlight here the importance, for this study, of the approximation we have made to construct a realistic atmosphere around a magnetic configuration that contains flux ropes (given by Eqs. (2), (3) and (7)). This allows us to account for the diversity of density in the coronal environment taking into account the importance of the magnetic field in shaping the structures in the solar corona, including dense horizontal structures.

An essential result from our modelling work is that the propagation of the ejection and the consequent perturbation at 4 R_⊙ is highly sensitive to the local pre-eruption conditions as the two analysed ejections carry different kinematic and plasma properties. Both also show no simple connection between the flux rope configuration at the lower boundary and its following evolution. Due to this, each flux rope ejection should be considered a unique phenomenon and consequently its associated space weather perturbation will also be unique.

In the present paper we have shown the potential capability of this technique through a specific example. An essential follow up to this work will be to compare models of magnetic flux rope ejections with observations. The quasi-static non-potential global model has been thoroughly tested for the identification of actual magnetic flux rope formation and times of eruptions. In the future we will also compare the magnetic flux rope evolution as described by the MHD simulation with observations of specific events.

5.2 Producing an operational model

The present model shows many interesting features that with further development may be a useful tool for application to space weather prediction. The technique presented could provide realistic time dependent inhomogeneous boundary conditions that couple the outer corona to the heliosphere. This coupling still requires development and testing within a heliospheric model. In particular, we have shown that from a single magnetic configuration we may produce multiple eruptions, where each of these eruptions have an inhomogeneous structure to the magnetic field and density as they pass through the interface between the corona and interplanetary space. Another interesting outcome is that the magnetic field and density properties of the two eruptions have significant differences in size and distribution. Therefore to some degree each eruption is unique, based on the initial properties of the flux rope and the surrounding non-potential field produced in the quasi-static evolution model.

For practical applications any space weather forecasting tool needs to have two key qualities: accuracy and efficiency. Accuracy is required in order to avoid false warnings and to ensure key parameters such as the arrival time of the CME at Earth and the orientation of this magnetic field and density are reliably determined. In terms of this our model sets high standards of accuracy as firstly it evolves the global corona from magnetic configurations determined from the time evolution of observed magnetograms. This produces non-potential magnetic configurations very close to observations and the global non-potential coronal model has proved highly reliable in predicting the locations, helicity and timing of flux rope formations, especially after the model has taken into account the emergence of new magnetic flux. Secondly, it solves the full MHD equations to follow the dynamical eruption itself. For time critical predictions efficiency is essential to run useful simulations that can output results before or close to when an event has taken place. In our tests, the global non-potential model is more than 10⁴ times computationally faster than a full MHD simulation. This is mostly due to the lower number of equations to solve and to the higher time steps involved when Alfvén waves do not need to be resolved. If a single MHD model was used, then the vast majority of the physical time of the simulation would be during the formation of the flux rope (days to weeks) compared to hours for the eruption. Thus the use of the global non-potential model to produce the initial eruptive configuration is twofold, first it solves fewer equations and secondly it is much more efficient. This means that only a small part of the physical time (once the ejection sets in) needs to be considered and modeled with a full MHD simulation. The combined model is in principle capable of predicting flux rope formation and subsequent onset of ejection days before it occurs, where the MHD simulation would only need to simulate a few hours of physical time which is well within current computational capacity. While this is a positive aspect of this approach, before such a technique can become viable as a semi-autonomous, factually accurate, operational space weather forecasting model, a number of practical and scientific obstacles need to be overcome. These include

• Flux emergence accuracy. The current global non-potential model is highly dependent on a semi-automated determination of flux emergence from observations. This is key as the emergence of new flux can significantly alter the magnetic configuration of the solar corona. Therefore a more accurate treatment must be developed. Techniques that allow for the direct assimilation of normal component magnetograms across all longitudes and latitudes are under development and testing (see Weinzierl et al. 2016a and Weinzierl et al. 2016b). One issue with the direct assimilation of magnetogram data into the simulation, is that only a portion of the global photospheric field may be observed at any given time. However, if magnetogram data were also available from an L5 based magnetograph this would significantly enhance the accuracy of the models. In the recent paper by Mackay et al. (2016) it has been shown that such L5 observations could increase the accuracy of global quantities relevant to eruptions on the Sun by anywhere from 26–40%. Pevtsov et al. (2016) reports similar improvements for the estimations of solar wind and plasma outflows as well.

• Identification of erupting flux ropes. Although the global non-potential model is able to describe the formation of magnetic flux ropes, it cannot account for the whole physics of an eruption. For the latter MHD simulations are employed. Therefore during the global non-potential simulations automated quantitative analysis of the properties of the coronal magnetic field must be carried out to first identify the existence of flux ropes and then the onset of ejections. A number of candidates for this task include instability analysis, measurement of thresholds, or recognizing critical indices in the parameters of the magnetic field configuration.

• Coupling AMRVAC back to the Global Non-potential Model. One of the main strengths of the global non-potential modelling technique is the continuous nature of the simulations over days to years, in which the transport of magnetic flux and magnetic helicity from low to high latitudes is followed. Once magnetic field configuration are imported from the global non-potential model into ARMVAC the continunity of the global non-potential simulations ends. Thus in order to allow a continuous modelling of the solar coronal evolution after the eruption, it is vital to couple the related MHD solution back to the global non-potential model. This will then be used to start a new segment of the quasi-static evolution. Currently this is not possible due to technical limitations where the primary variable in the Global Non-potential Model is the vector potential A, while in MPI-AMRVAC it is the magnetic field B. To resolve this issue a new MPI-AMRVAC module that solves the MHD equations using A is under development.

• Matching the simulations with the lower boundary of space weather models. Most space weather forecast tools focus on the physical domain extending from 0.1 AU to > 1 AU. The lower boundary of this range corresponds to about 21 R_⊙ and is significantly higher than the outer boundary of the simulation shown here. To bridge the gap between the models a number of options will be investigated in the future. These include using a simple solar wind model above 4 R_⊙ (for a review: Echim et al., 2011) or to scale the plasma and velocity fluxes from 4 R_⊙ up to 21 R_⊙. In the latter approach we would consider it as a perturbation to a background solar wind that is included in the applied space weather forecast tool. In future studies we will consider the validity and technical aspects of both of these approaches.

In conclusion, the present work shows that in the future the coupling of the quasi-static non-potential evolution model with MHD simulations may become a useful and essential tool to bridge magnetogram observations with space weather modellling. This is due to the complexity of the magnetic configurations generated in the solar corona and the interaction between propagating flux ropes and the coronal environment. The technique we put forward is a viable one for future space weather forecasting tools and has the potential not only to improve our forecasting capabilities, but also to lead to a better understanding of the physics of CME origin.

Acknowledgments

This research has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 647214). This work used the DiRAC Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk. This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, STFC capital grant ST/K00087X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. We acknowledge the use of the open source (gitorious.org/amrvac) MPI-AMRVAC software, relying on coding efforts from C. Xia, O. Porth, R. Keppens. D.H.M. would like to thank STFC and the Leverhulme Trust for their financial support. ARY was supported by STFC consortium grant ST/N000781/1 to the universities of Dundee and Durham. The editor thanks two anonymous referees for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 Representation of the magnetic configuration imported from the quasi-static non-potential evolution model. The two panels show maps of the radial component of the magnetic field on the surface from two opposite points of view. Blue lines are some representative magnetic field lines to highlight existing structures.
In the text

	Fig. 2 Mollweide projections (central meridian at longitude ϕ = 180∘ ) on the lower boundary of the MHD simulation of (a) Lorentz Force, (b) ω, (c) Log10(ρ[g∕cm3]) and (d) Log10(T[K]).
In the text

	Fig. 3 Maps of the column density seen from either side of the Sun in the MHD simulation along with some superimposed magnetic field lines traced from the lower boundary (green) at t = 0 min (a), t = 33.6 min (b), and t = 60.3 min (c). Central meridian for left-hand side panels is at ϕ = 90∘ and at ϕ = 270° for right-hand sided panels and the FOVs are chosen to be same as in Figure 1.
In the text

	Fig. 4 (a) Cuts of Log10(ρ[g∕cm3]) along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line). (b) Cuts of Log10(T[K]) along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line).
In the text

	Fig. 5 (a) Cuts of radial velocity along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line). (b) Cuts of Log10(ω) along the radial direction from solar surface to the outer boundary at the locations of FR1 and FR2 at t = 10.4 min (blue line) and t = 27.8 min (red line).
In the text

	Fig. 6 Energy in the MHD simulation integrated over the whole spatial domain as a function of time. Black line shows total energy, blue magnetic energy, red kinetic energy, and green thermal energy.
In the text

	Fig. 7 Mollweide projections (central meridian at longitude ϕ = 180°) at the outer boundary of the MHD simulation at t = 60.3 min of (a) ρ[g∕cm3] and (b) radial velocity of the plasma [km∕s].
In the text

	Fig. 8 Cuts of (a) Log10(ρ[g∕cm3]), (b) radial velocity [km/s], and (c) temperature [K] along ϕ at t = 0 min (black line), t = 33.6 min (blue line) and t = 60.3 min (red line) atthree different co-latitudes: the center of FR1 (θ = 71.4°), the equator (θ = 90°), and the center of FR2 (θ = 110°).
In the text

	Fig. 9 Mollweide projections (central meridian at longitude ϕ = 180°) at the lower boundary of the MHD simulation at t = 0 min (left column) and at the outer boundary at t = 0 min (central column) and t = 60.3 min (right column) for the three magnetic field components Br, Bθ, Bϕ.
In the text

	Fig. 10 Cuts of Bθ along ϕ at three different latitudes: the center of FR1 (θ = 71.4°), the equator (θ = 90°), and the center of FR2 (θ = 110°). a) At t = 0 and lower boundary. b) At t = 0 min (black line), t = 33.6 min (blue line) and t = 60.3 min (red line) and outer boundary
In the text

	Fig. 11 Mollweide projection of ω at the outer boundary of the MHD simulation at t = 60.3 min (central meridian at longitude ϕ = 180°) .
In the text