© J. Pomoell and S. Poedts, Published by EDP Sciences 2018

1 Introduction

From a socio-economic perspective, among the most important manifestations of the dynamic nature of the Sun are solar eruptions. In particular coronal mass ejections (CMEs) − large magnetized clouds of plasma that erupt frequently from the Sun − constitute the main drivers of adverse conditions in the inner heliosphere (e.g., Richardson et al., 2001). The detectable space weather effects on Earth appear in a broad spectrum of time and length scales and have various harmful effects for human health and for our technologies on which we are ever more dependent (e.g., Hapgood, 2011; Green & Baker, 2015; Schrijver et al., 2015; and references therein).

Severe conditions in space can hinder or damage satellite operations as well as communication and navigation systems and can even cause power grid outages leading to a variety of socio-economic losses. Following a growing awareness of these effects, several studies have attempted to estimate the cost of a severe space storm. A report from the National Research Council of the National Academies in 2008 (Press, 2008) estimates that the advent of an event such as the one which occurred in September 1859 (the Carrington event) would cost today on the order of one trillion (10¹²) USD and would take 4–10 years in repairs to recover, an order of magnitude more than that of Hurricane Katrina (Allianz, 2009). More recently, Eastwood et al. (2017) estimated that the total economic loss varies between 0.5 and 2.7 trillion USD based on calculations examining disruption to the global supply chain. With an alternative methodology they found a total loss of 140–613 billion (10⁹) USD.

The average occurrence of extreme events at the level of the Carrington event is estimated to one in 250 years with a confidence level of 95% (Cannon et al., 2013). Riley (2012), however, estimates that there is a 12% probability of such a super storm occurring within the next decade. The “normal” i.e., noncatastrophic, space weather events, however, also have a considerable impact. Information from insurance claims, for instance, suggests that geomagnetically induced currents cause losses to the U.S. power grid alone may be of the order of 5–10 billion USD per year (Schrijver et al., 2014). Hence, NASA, ESA and other space agencies are seeking to develop models of space weather events and to develop tools to reliably predict effects on Earth so that protective measures can be taken.

Despite decades of research, the accuracy of long-term space weather forecasting is still at best modest. A significant focus of the community has been on predicting the arrival time of CMEs. Reviewing a number of methods, both empirical and physics-based, Zhao & Dryer (2014) find that the root-mean-square error in the prediction of the arrival time generally is 12 h and 10 h for the mean absolute error. Similar errors are reported by Mays et al. (2015) who utilize ensemble modeling employing the WSA–ENLIL model.

While predicting the time of arrival is certainly crucial, a major complication is that not all CMEs are equally geo-effective, and thus, it is paramount to improve our ability to predict their geo-effectiveness, see e.g., Siscoe (2007); Zheng (2013); Lavraud & Rouillard (2013). The main problem is that our ability to predict the magnetic field characteristics of CMEs, a major factor determining their geo-effectiveness, is minimal. Only recently have space weather oriented methods that aim to predict the magnetic structure of CMEs been constructed (e.g., Savani et al., 2015; Isavnin, 2016; Kay et al., 2017). A significant hurdle for such empirical-based methods is to account for the dynamics of the CME propagation especially for complex cases that include interacting solar wind structures. One viable approach is to employ physics-based modeling, e.g., magnetohydrodynamic simulations, that self-consistently capture the complex dynamics. At the moment, the ENLIL model is the only MHD heliospheric model that is operational (Parsons et al., 2011), i.e., available and routinely used for forecasting.

In this paper, we describe the current implementation of the European heliospheric forecasting information asset (EUHFORIA). This newely developed model is a physics-based simulation tool designed for space weather forecasting purposes. The modeling described in this paper employs proven semi-empirical models and methods to compute the boundary conditions that drive the magnetohydrodynamics model of the inner heliosphere that forms the core of EUHFORIA. While this paper focuses on detailing the baseline implementation, features that reach beyond the baseline will be presented in future work. Furthermore, a more detailed validation, considering several historic cases and a wider simulation parameter range, will also be conducted later.

This paper is structured as follows. Section 2 contains the technical description of the different parts of the model. Section 3 contains an initial validation test case. In Section 4 a discussion of the results and an outlook on upcoming future improvements of the model are considered.

2 EUHFORIA modeling scheme

The spatial domain of EUHFORIA is divided into two distinct regions: a coronal domain extending from the photosphere or low corona up to a heliocentric distance of r = R_b = 0.1 AU, and an inner heliosphere domain from 0.1 AU typically extending up to 2 AU so as to keep the orbit of Mars well inside the domain. The premise of the division is that at 0.1 AU and beyond the solar wind plasma is supersonic and super-Alfvénic which means that no information is traveling towards the Sun (all MHD characteristic curves are outgoing). As a result, the heliosphere model is agnostic of the coronal model, and makes coupling the two domains a one-way process. This allows using any available coronal model for the coronal domain, independent of whether the model is empirical or physics-based. Indeed, this property has been employed for instance in the context of the heliospheric ENLIL model that accepts data from different coronal models (Owens et al. (2008); Gressl et al. (2013)). Also MHD simulations of the entire corona-inner heliosphere system frequently employ a division of the two domains, see e.g., Lionello et al. (2013) and references therein.

2.1 The coronal domain: semi-empirical solar wind model

The essential task of the coronal model in EUHFORIA is to produce plasma conditions at r = 0.1 AU that provide the boundary conditions needed by the inner heliosphere model. The model should provide realistic values for the MHD variables that represent the prevailing large-scale solar wind streams as well as transient disturbances due to propagating CMEs. In the current modeling scheme these two tasks are separated. The method by which CMEs are included is described in Section 2.3, while the undisturbed ambient solar wind model is detailed next.

The solar wind model of EUHFORIA relies on semi-empirical relationships between topological properties of the coronal magnetic field and the measured solar wind parameters. Such an approach was selected in favor of a more complete physics-based approach such as a coronal MHD model by considering the accuracy of the model in relation to the computational requirements. Comparative studies have shown that empirical models employed as means to provide boundary conditions for heliospheric models perform well when compared with physics-based coronal models (e.g., Owens et al. (2008); Gressl et al. (2013); Jian et al. (2015, 2016)).

The empirical model employed in this work follows the path laid out by the Wang–Sheeley model (Wang & Sheeley, 1990), the DCHB model (Riley et al., 2001) as well as the Wang–Sheeley–Arge model (Arge, 2003) in that the solar wind speed is given as a function of the flux tube expansion factor f and the distance of the foot-point of the flux tube to the nearest coronal hole boundary d. As such, the wind speed is determined completely by the properties of the global three-dimensional coronal magnetic field.

2.1.1 Magnetic field model

Similar to Arge (2003); Detman et al. (2006), the magnetic field model consists of two parts: in the lower corona for r ∈ [R_⊙, R_i], the magnetic field is given by the potential field source surface (PFSS) model (Altschuler & Newkirk, 1969) in which the magnetic field is assumed to be current-free and is set to be purely radial at a given source surface radius R_ss ≥ R_i.

In the upper coronal domain, r ∈ [R_i, R_b], the magnetic field is given by the Schatten current sheet (SCS) model (Schatten et al., 1969), which involves solving an additional Laplace equation similar as in the PFSS model but with different boundary conditions. The absolute value of the radial field given by the PFSS model at r = R_i, i.e., |B_r;PFSS(R_i, θ, φ)|, specifies the input data for a Neumann boundary condition for the magnetic scalar potential at the interface between the PFSS and SCS model, while at infinity the SCS field is required to vanish. The main purpose of the SCS model is to extend the magnetic field in a nearly radial fashion while retaining a thin structure for the heliospheric current sheet. Note that if the source surface radius is set to be the same as the interface radius, R_ss = R_i, the magnetic field lines do not make a smooth transition at the boundary of the two models since the boundary conditions are not compatible. To minimize this effect while at the same time approximately preserving the amount of flux open to the solar wind, we set R_i = 2.3R_⊙ and R_ss = 2.6R_⊙ as the default values as suggested in McGregor et al. (2008). The Laplace equation appearing in both the PFSS and SCS models are computed using an expansion in solid harmonics. By default, the expansion is computed up to degree l = 120 (l = 70) in the PFSS (SCS) model, although the parameters can be set at run-time by the modeler.

With the three-dimensional magnetic field determined, the large-scale topology of the coronal magnetic field is characterized by tracing a set of field lines through the model. Specifically, to locate the regions in the low corona that are open to the solar wind (i.e., coronal holes), a field line starting from the base of the model is traced upwards (away from the Sun) for each pixel on a map that covers the spherical surface with a given angular resolution. If the point at which the tracing stops is high in the corona, the pixel is designated open, while if the field line returns to the corona the pixel is designated as belonging to a region with closed magnetic topology. To compute the connectivity of open field lines to their coronal source regions, a second set of field lines are traced down to the corona for each pixel on a spherical surface located at the outer boundary (0.1 AU) of the model. This allows to compute the flux tube areal expansion factor f: (1) where B_r(R_p, θ_p, φ_p) is the radial magnetic field at a point p in the upper corona at radius R_p = R_b and B_r(R_⊙, θ, φ) the radial magnetic field at the photospheric footpoint of the field line passing through p. Note that we define the flux tube expansion factor using only the radial magnetic field component following Riley et al. (2015) and is different from other works (e.g., Arge & Pizzo (2000)) that instead use the magnitude of the magnetic field vector. While at the outer boundary the magnetic field in the SCS model is predominantly radial, at the photospheric boundary the horizontal components are significant.

In addition to the flux tube expansion factor, the great-circle angular distance d from the foot point of each open field line to the nearest coronal hole boundary is computed.

2.1.2 Empirical solar wind plasma specification

The task of the empirical solar wind model in EUHFORIA is to provide the plasma parameters characterizing the ambient solar wind at 0.1 AU for a given flux tube. The quantity around which the empirical model is built is the solar wind speed v, and is obtained through a given function of the magnetic field properties: v = v(f, d). Several different forms for specifying the wind speed have been presented in the literature (e.g., Arge & Pizzo (2000); Arge (2003); Riley et al. (2001); Owens et al. (2008); Detman et al. (2006); Wiengarten et al. (2014) McGregor et al. (2011)).

In EUHFORIA, the functional form is up to the choice of the modeler, and is specified as an input parameter when running the model. In this work, we employ the following prescription: (2) with parameters v₀ = 240 km/s, v₁ = 675 km/s, α = 0.222, β = 1.25 and w = 0.02 rad. The functional form as well as the parameter values, except for w, are identical to that given in van der Holst et al. (2010) and McGregor et al. (2011) (the “original WSA model” as given in their equation (2)). The parameter w serves to normalize the dependence of the solar wind speed on the distance to the nearest coronal hole boundary.

In the original WSA approach, the empirical formula such as equation (2) specifies the solar wind speed at Earth since the solar wind is propagated to Earth assuming a constant velocity but with stream interactions taken into account. In the MHD approach, on the other hand, the wind speed at 0.1 AU is required instead. Since the solar wind continues to accelerate beyond 0.1 AU in the MHD model, we subtract from equation (2) a constant value of 50 km/s in order to avoid systematically overestimating the wind speed. Furthermore, the resulting speed is capped to be in the range v_r ∈ [275, 625] km/s (see McGregor et al. (2011)). Finally, the obtained solar wind speed map is rotated by 10^∘ to account for solar rotation that is not included in the magnetic field model.

The empirically determined speed of the ambient solar wind forms the basis for prescribing the boundary conditions at 0.1 AU for the heliospheric MHD model. As the plasma is super-fast, all eight MHD variables need to be specified. Following, the velocity is chosen to be purely radial, with the radial component v_r set equal to the empirical speed prescription. For the magnetic field, the meridional component is set to zero. Then, for a given radial component B_r, the azimuthal component is specified as B_φ = − (B_r/v_r)R_bΩsinθ resulting in a zero electric field in the co-rotating frame which is consistent with the requirement to obtain a steady-state solution in the co-rotating frame. The radial component itself is computed as a function of the wind speed, (3) where sgn(B_corona) is the sign of the magnetic field as given by the coronal model, v_fsw = 675 km/s refers to the speed of the fast solar wind that carries a magnetic field of B_fsw = 300 nT. We opt to determine the magnitude of the magnetic field not from the coronal magnetic field model but rather from the computed wind speed in order to avoid the so-called open flux problem, i.e., the issue that the magnetic field strength in interplanetary space is often underestimated when inferred from coronal models, see e.g., Linker et al. (2017) and references therein.

The plasma number density is given by (4) with n_fsw = 300 cm⁻³ the number density of the fast solar wind. This prescription ensures a constant kinetic energy density on the spherical surface at r = R_b. Finally, the plasma thermal pressure is chosen to be constant on the boundary and equal to P = 3.3 nPa which corresponds to a temperature of T_fsw = 0.8 MK in the fast solar wind. Whereas similar choices for the plasma density and pressure have been employed (e.g., Odstrcil & Pizzo (1999)), other options have also been explored. For instance, Detman et al. (2006) set the density inversely proportional to v_r thereby enforcing a constant mass flux and further assume the total pressure (thermal plus magnetic) to be constant. In another approach, Shiota et al. (2014) set the density and temperature utilizing empirical fits to Helios in situ data as computed by Hayashi (2003).

2.1.3 Magnetogram input

The fundamental component of the semi-empirical coronal model described above is the magnetic field model. For providing the synoptic magnetograms required by the PFSS model, data provided by the Global Oscillation Network Group (GONG) is employed. In particular, the hourly updated standard synoptic magnetograms are used by default. The maps are projected to Plate Carée with a 1° resolution and smoothed using a Gaussian filter with standard deviation equal to 0.8. In Figure 1 the process of constructing the semi-empirical coronal model is summarized in the form of a flowchart.

Fig. 1 Schematic of the steps involved in constructing the semi-empirical coronal model that provides the boundary conditions for the heliospheric MHD model.

2.2 The inner heliosphere domain: MHD model

From R_b = 0.1 AU onwards, the inner heliosphere model in EUHFORIA consists of a three-dimensional time-dependent magnetohydrodynamics simulation that self-consistently models the propagation, evolution and interaction of solar wind streams and CMEs. The equations solved are those of ideal MHD augmented with gravity: (5) (6) (7) (8) where ρ, v, B, E, P are the mass density, velocity, magnetic field, total energy density and thermal pressure of the plasma, respectively. The total energy density is given by where γ is the polytropic index, while the gravitational acceleration with G the gravitational constant and M_⊙ the solar mass. The equations are solved in a frame that corresponds to the Heliocentric Earth Equatorial (HEEQ) coordinate system. As such, the inner radial boundary data rotates with respect to the computational grid at the solar rotation rate. Although the chosen frame is not inertial, we omit the Coriolis and centrifugal terms arising as a consequence of the orbital motion of Earth as their contribution is small.

The polytropic index is chosen equal to 1.5 as in Odstrcil et al. (2004). The reduced index is a simple way of modeling additional heating and thereby acceleration of the solar wind (see, e.g., Pomoell & Vainio (2012) for further details).

The numerical scheme to solve the MHD equations is a finite volume method together with a constrained transport approach for advancing the magnetic field in an exactly (up to machine accuracy) divergence-free way. In the method, the hydrodynamical quantities (ρ, v, E and P) are cell volume-averaged quantities, whereas the magnetic field components are cell face-averaged and the electric field components are cell edge line-averages and are co-located in a staggered fashion on the grid. Standard piece-wise linear reconstruction and an approximate Riemann solver are employed to arrive at a robust second-order accurate scheme (see Kissmann & Pomoell (2012); Pomoell & Vainio (2012)). The computational domain extends from 0.1 to 2 AU in the radial direction and spans 120° in latitude and 360° in longitude. The mesh used in this work is uniform in all directions with a 4° angular resolution and a radial grid spacing of Δr ≈ 0.0074 AU ≈1.6R_⊙ (256 cells in radial direction).

To implement the empirically determined solar wind boundary conditions in the MHD model, two layers of ghost cells are used the values of which are determined at each time step through a linear extrapolation using the boundary condition data (with the boundary data at r = R_b located at the interface between cells) and the data in the first in-domain cell. Before this procedure, the rigidly rotating boundary map is interpolated to the cell faces at the boundary interface. A complication lies in the fact that the constrained transport method requires the tangential electric field components to be specified at the boundary. Similar to Lionello et al. (2013), the tangential electric field is determined using a poloidal-toroidal decomposition (PTD) by writing (9) where ψ is a scalar potential. The evolution of the radial component of the magnetic field on the spherical surface is then given by (10)

Thus, for a given time-evolution of B_r, the electric field can be computed by solving the Poission equation for the scalar potential ψ. Since we assume the inner boundary to rotate rigidly, the left hand side is simply a rotation in longitude. In fact, in this simple case, the tangential electric field can be computed analytically as . Despite this simplification in the current work, we utilize the PTD machinery with future applications in mind. Note that in general, a gradient term should also be present in equation (9). We set it to zero in the co-rotating frame which is consistent with the chosen boundary conditions of the magnetic field (see Sect. 2.1.2). At the outer radial boundary, open boundary conditions implemented via simple extrapolation is used, whereas at the latitudinal boundaries symmetric reflection is employed.

2.3 The CME model

CMEs constitute significant global scale transient structures in the solar wind. To model their influence on the plasma in the inner heliosphere, EUHFORIA employs the cone CME model, similar to that of Odstrčil & Pizzo (1999). The model treats the CME as a hydrodynamic cloud exhibiting a simple geometry as the CME evolves in the upper corona characterized by a constant angular width, propagation direction and speed. The cross section is assumed to be circular, and the CME is filled homogeneously (i.e., density, pressure and radial speed are constant).

The CMEs are introduced as a time-dependent boundary condition at the inner radial boundary at 0.1 AU. Following the assumption of a circular cross-section, the perturbation on the boundary has a circular shape with a radius that varies in time. The speed at which the radius of the circular enhancement increases and decreases is linked to the speed and the shape of the CME. In this work, we set the angular width of the enhancement at 0.1 AU to be sinusoidal in time which corresponds approximately to pushing a sphere through a plane at constant speed. Specifically, the ambient solar wind values for each quantity q ∈ {ρ, v_r, T} are replaced by constant values {ρ_CME, v_CME, T_CME} for points at r = R_b for which (11) where the width as a function of time is (12)

In the equations above, θ_CME, φ_CME are the co-latitude and longitude of the propagation direction of the CME center, t_onset is the time at which the CME reaches r = R_b and t_half = R_b tan(ω_CME/2)/v_CME with ω_CME the angular width of the CME.

The kinematic properties of a given CME can be estimated using coronagraph data and fitting tools such as StereoCat (Mays et al., 2015) or the forward modeling tool of Thernisien et al. (2009). In addition to the kinematic parameters, the density and temperature of the plasma cloud need to be estimated. Thus, the parameters that the forecaster is required to supply for each CME are the mass density (ρ_CME), temperature (T_CME), velocity (v_CME, θ_CME, φ_CME), angular width (ω_CME) and onset time at 0.1 AU.

2.4 Phases of a run

A run of the heliospheric MHD model consists of thee distinct successive phases, as depicted in the schematic of Figure 2. Each simulation starts by performing a relaxation in which the MHD equations are advanced in time starting from a suitable initial condition. The purpose of this phase is to obtain a steady-state (in the co-rotating frame) MHD solar wind solution that is consistent with the provided boundary conditions and that forms the proper starting point of the model. As a result, the simulation is not sensitive to the initial guess for the heliospheric plasma state. By default, a relaxation time of 14 days is employed which is sufficient for a wind with a speed of ∼250 km/s to traverse 2 AU and exit the computational domain through the outer radial boundary.

At any given time, at least several CMEs are present in the inner heliosphere on average, as even during solar minimum several CMEs are launched per day (Robbrecht et al., 2009). The plasma environment in the interplanetary space is significantly modified by the presence of the CMEs, and this needs to be taken into account when modeling the evolution of a particular eruption. To accomplish this, after the solar wind solution has been obtained, the model enters a phase where CMEs that have occurred in the previous days are inserted. Typically, significant CMEs from the five days prior to the start of the forecast are inserted.

The last phase of the run is the actual forecast. In EUHFORIA, the starting time of the forecast is defined equal to the observation time of the magnetogram that is used to construct the coronal model. The forecast is reasonable to run forward in time approximately up until the East limb observation crosses the central meridian, and corresponds to ∼5–7 days.

Fig. 2 The phases of the heliosphere simulation with typical durations. Each run starts with a solar wind relaxation, followed by a phase in which previously observed CMEs are introduced. The forecast starts at t = 0 and is typically run for 5 days.

3 Initial validation run: the events of June 17–29, 2015

As a first validation test, in this section, we present results of applying EUHFORIA to model the plasma conditions in the inner heliosphere during June 17–29, 2015. This particular time period was chosen as it was the first event for which a forecast was produced using an early version of the model. To conduct the validation run, the modeling pipeline is used completely similarily as when run in operational mode, i.e., the boundary conditions (solar wind and CMEs) are specified identically as would be done when producing a forecast. In particular, the in situ data that in a true operational setting would not be available, is in no way used to optimize the boundary conditions.

3.1 Overview of events

A succession of solar eruptions occurred during the time period from June 17 to June 29, 2015. Most of the on-disk activity was caused by NOAA AR 12371 which produced three GOES M-class flares and numerous C-class flares. The on-set time of the strongest flare (M7.9) occurred 08:02 UT on June 25. A number of CMEs were also observed. The CACTUS catalog, which uses automated detection of transients in SOHO/LASCO coronagraph imaging observations Robbrecht et al. (2009), lists 51 events for the time period. Seven distinct events are marked as potential (partial or full) halo CMEs, with the three that have a plane-of-sky angular width larger than 270 closely related in time with the M-class flares from AR 12371.

In the vicinity of Earth, six interplanetary shocks for the time window are cataloged in the Heliospheric Shock Database (www.ipshocks.fi) Kilpua et al. (2015), all observed by the Wind spacecraft. The events were coincident with a geomagnetic storm that reached a peak Dst of −204 nT on 05 UT, June 23 as given by the provisional Dst index data provided by the World Data Center of Geomagnetism in Kyoto, Japan.

3.2 Coronal model

To compute the coronal model, the input magnetogram must first be selected. As noted at the beginning of this section, we run EUHFORIA in a simulated operational mode by considering a scenario where the CME that occurred on June 25, 2015 (CACTUS on-set time 08:36 UT) has just been observed and a forecast of the impact of the eruption on the heliospheric conditions is to be computed. Thus, a GONG magnetogram prior to the eruption (01:04 UT) is selected as the input to the coronal model.

Figure 3 presents the radial magnetic field component B_r at the inner radial boundary of the coronal magnetic field model (PFSS). Note that the plotted data is not the raw input magnetogram, but rather the radial magnetic field component that is given by the PFSS model after the magnetogram has been projected and processed (see Sect. 2.1.3). As a result, Gibbs-type artifacts due to the finite number of harmonics are noticeable especially around the strongest active region fields. The large scale structure of the magnetogram is, however, well reproduced. In Figure 4, the open and closed magnetic field regions as determined by tracing field lines (Sect. 2.1.1) are shown. In the image, gray areas correspond to closed field regions (lighter grays indicating fields pointing away from the Sun at the footpoints) while blue (red) depicts regions of open field lines with the magnetic field directed towards (away from) the Sun.

Fig. 3 Radial component of the magnetic field as given by the coronal magnetic field (PFSS) model at r = R⊙.

Fig. 4 Open and closed regions as determined by the coronal magnetic field model. Blue (red) pixels correspond to open field lines that continue in to the solar wind with a magnetic field pointing toward (away from) the Sun. Gray pixels indicate regions of closed magnetic topology, with lighter gray indicating a positive field at the foot point.

3.3 Heliosphere model

While the empirical coronal model provides the state of the ambient solar wind, CMEs need to be inserted separately in order to model the disturbed heliospheric conditions during times of significant eruptions, as is evidently the case for the chosen time interval. As discussed in Section 2.3, an accurate determination of the cone model parameters requires fitting the geometric model to a sequence of coronagraph observations.

In order to maintain objectivity in the results of the modeling, we use the cone model parameters provided by the CCMC Space Weather Database Of Notifications, Knowledge, Information (DONKI) database directly. For the given time period the DONKI database lists in total 26 CMEs, 21 of which are slow (less than 500 km/s) or narrow (half-width less than 35 degrees) or not directed towards Earth (source longitude not within ±60 deg HEEQ). The remaining 5 CMEs are selected for modeling. The selected CMEs exhibit a clear halo signature, and events 1, 4 and 5 are related to the M-class flares from AR 12371. The parameters of the five selected CMEs are provided in Table 1.

While the DONKI database entries list the kinematic parameters of the cone model fittings, the density and temperature of the modeled CMEs are not provided. For simplicity, we assume each CME to have the same density ρ_CME = 10⁻¹⁸ kg m⁻³ and temperature T_CME = 0.8 MK.

3.3.1 Heliospheric dynamics

Figures 5–8 present four snapshots from the MHD simulation depicting the heliospheric dynamics during the events. Movies of the dynamics is available in the electronic journal. In each figure, the top two panels represent the radial speed in the simulation (in units of km/s), while the lower two panels plot the number density (in units of cm⁻³) scaled by r² measured in astronomical units, (13)

In each figure, the left panel plots the quantity in the heliographic equatorial plane, while the panels to the right show the values in a meridional plane including Earth. The positions of the inner planets as well as the locations of the twin STEREO spacecraft are indicated with the markers. Due to the HEEQ coordinate system, the location of Earth remains fixed at the same longitude, but the latitude and radial distance vary in time. The gray straight lines are drawn as guides for the eye, while the circles are drawn at heliocentric radii of r = 0.5, 1, 1.5, 2 AU.

The figures reveal the near-Earth environment to be highly dynamic as a result of the successive eruptions. While the CMEs drive shocks ahead of them, significant low-density regions are formed in the wakes of the CMEs. As a result, the CMEs (except for CME1) are propagating in disturbed conditions caused by the preceding eruptions.

As is seen in Figure 5, CME1 is launched in a direction significantly to the East of the central meridian. Only the flank of the CME reaches the position of Earth. This is also evident in the meridional slice of Figure 5 where a relatively weak signal is observed. The propagation of CME1 reveals that the ambient solar wind structure impacts the evolution of the CME considerably. CME1 interacts both with slow and faster wind streams. In the region of faster wind and lower ambient density, the plasma of the CME is compressed considerably less than at the parts of the front propagating through slower wind. This causes the CME front to become dimpled, with a local minimum of density as indicated by the black arrow in Figure 5.

In contrast to CME1, CMEs 2–4 are all on a direct trajectory towards Earth. Although the launch of CMEs 2 and 3 are separated in time by 38 h, they interact considerably as CME3 is approximately twice as fast as CME2. According to the model, they reach Earth at nearly the same time, with CME2 preceding CME3 by only a few hours, as is seen in Figure 6. The dynamics at this time (around 11 UT on June 22) is especially complex as there are four structures that reach the position of Earth within a time window of roughly 24 h: a high-density solar wind stream, the flank of CME1, the slow CME2 and the overtaking faster CME3. Soon after passing Earth, CMEs 2 and 3 merge and propagate as a single structure. CME4, launched 40 h after CME3, propagates in the plasma wake of CMEs 2 and 3. In spite of the relatively homogenous upstream plasma conditions that do not show any clear solar wind streams, the front of CME4 becomes non-uniform beyond 0.5 AU when it reaches the density minimum of the wake of the preceding CMEs. In Figure 7, especially the density front shows a minimum at the leading edge of the CME that roughly coincides with the Sun–Earth line.

CME5, inserted at 0.1 AU at 10:51 UT on June 25, is clearly launched to the west of the central meridian (φ = 46^∘ HEEQ). Until roughly 03 UT on June 27 when the nose of the CME has reached ∼1 AU, the CME appears to be on a trajectory that would not intercept Earth. However, at this time the eastern flank of the CME expands rapidly in the lateral direction, coincident with the front reaching a region of low density (seen in Fig. 8 as the purple region surrounding the near-Earth space) of the wake of CME4. This lateral expansion causes the CME to produce a signal at Earth.

Fig. 5 Snapshot of the MHD simulation at 03:03 UT on June 21, 2015. Top row shows the radial speed, while the bottom row shows the scaled number density. The left panels depict the solution in the heliographic equatorial plane, while the right panels show the meridional plane that includes Earth.

Fig. 6 Same as Figure 5 but at time 11:01 UT on June 22, 2015.

Fig. 7 Same as Figure 5 but at time 11:01 UT on June 24, 2015.

Fig. 8 Same as Figure 5 but at time 11:03 UT on June 27, 2015.

3.3.2 Comparison with in situ observations

Figure 9 shows the radial speed at Earth as a function of time obtained from the simulation (with a cadence of 10 min) together with the speed from five-minute OMNI data. In addition, the time of the six shocks listed in the ipshocks database are shown as vertical gray dotted lines. Similarly as in Figure 9, in Figure 10, the number density of the in situ observations and simulation data is shown.

At the beginning of the simulated time period, the solar wind at Earth is characterized by a steadily decresing solar wind speed as the preceding high speed stream transitions to a stream of slow wind. As is evident in Figure 9, the model captures the decreasing solar wind speed remarkably well. Interestingly, the observed density of the wind does not show a clear transition, instead remaining at roughly 2–3 cm⁻³. This behaviour is contrary to the simulation results in which the density increases steadily as the wind speed decreases. From the beginning of the simulated period to ∼June 20, the simulated density follows closely a behaviour consistent with equation (4) after which the density increases more rapidly due to the interaction with the slow wind stream. A change in the behaviour in the density is visible in the observed data at ∼18 UT, June 20 at which time the density starts to increase.

At 16 UT on June 21, a fast-forward shock is observed in the WIND data (marked with S1). Approximately two hours later the simulation data shows a gradual increase of the plasma velocity, peaking at ∼01 UT on June 22. The increase is due to the flank of CME1 reaching Earth. The density data reveals the observed large-scale structure to be sharper than the simulated one, indicating the flank of the modeled CME to be too wide. At ∼06 AU on June 22, a second shock (S2) is observed and a third (S3) follows approximately 13 h later. In the simulation, two shocks are present later on June 22 that are separated only by 4 h. The first of these, indicated by the arrow in Figure 9, is the signal of the slow CME2 reaching Earth. Although the signature of CME2 appears in the simulated data at Earth ∼12 h later than S2, the similarity between the observed and simulated density and velocity jumps suggest that S2 is the shock driven by CME2. Thus, the simulated CME2 is delayed, either as a result of interaction with the dense solar wind (see Figs. 5 and 6) or a too low injection speed. On the other hand, the observed shock S3 and the shock driven by CME3 are nearly coincident, with their time of arrival within three hours. The density of the simulated CME is overestimated. However, it is difficult to determine the source of the error as disentangling the relative contributions of the interacting structures (see Sect. 3.3.1) is not possible without a study of the influence of e.g., the assumed density of the CME. It is interesting to note that although the simulated density following S3 is too large, the morphology is nevertheless remarkably similar.

Similar to CME3, the observed shock S4 and the shock driven by CME4 are very closely related in time (within ∼30 min of each other). This is the case even if the jump in velocity is overestimated in the simulation. On the other hand, the last simulated event CME5 produces a signature at Earth clearly late as compared to the increase in speed observed after the shock (S6). However, the simulation interestingly shows a small jump in speed and density coinciding in time with that of S6. Note that the density plot (Fig. 10) does not appear to show any increase after S6. At the time of S6, the OMNI data set features a data gap, whereas WIND data shows a jump in density. Nevertheless, the simulation clearly overestimates the density at this time. The structure responsible for the density bump is the high-density flank of CME5, indicated by the black arrow in the lower panel of Figure 8. However, the reason for this discrepancy between the simulation and observations is difficult to pinpoint. Closer to the Sun, the CME flank clearly interacts with a slow stream forming in the wake of CME4. It is possible that the wind stream forms too rapidly after CME4 (see also below). On the other hand, even the lowest-density region of the CME located at the leading edge along the ∼30° longitude line that does not encounter the slow wind is too high (∼10 cm⁻³). Thus, it is possible that the density in particular at the flanks of the CME are overestimated.

The behavior of the radial velocity in the wakes of CMEs 3 and 4 shows a similar behavior when contrasted with the observations: the speed decreases more rapidly towards the ambient solar wind speed in the simulation as compared with the observations. This would suggest that the wakes of the CMEs simulated by EUHFORIA are too narrow in radius. On the other hand, the density in the wake of CMEs 3 and 4 captures the decrease rather well. However, if the density of CME3 was not overestimated, it is likely that the density would drop below the observed values in the wake similar to the speed.

Fig. 9 Radial speed at the position of Earth as a function of time. The blue curve shows the data from the simulation (saved at 10 min cadence), while the red curve shows OMNI 5 min data.

Fig. 10 Number density at the position of Earth as a function of time. The blue curve shows the data from the simulation (saved at 10 min cadence), while the red curve shows OMNI 5 min data.

4 Conclusions and outlook

In this paper, we have presented EUHFORIA, a forecasting-capable physics-based simulation model of the inner heliosphere driven by boundary conditions based on empirical modeling methods. Following a detailed account of the simulation methodology in Section 2, an initial validation run focusing on the time period of June 17–29, 2015 was shown.

The validation run, featuring five cone-model CMEs, presented a highly dynamic heliosphere. This is especially the case at ∼11 UT on June 22 when three CMEs and a high-density slow wind stream interact and reach Earth within a short time span. In spite of the complex dynamics, the in situ observations from the near-Earth space and the simulation data agree well for this particular time-period, with the model predicting arrival times for two of the CMEs within three hours of the observed interplanetary shocks. However, it is important to note that drawing conclusions regarding the accuracy of the model for other events or time-periods cannot be made based on the sole event presented in this work. Comprehensive validation tests in which a large number of events are studied is required in order to quantify the accuracy of the model.

The results of the validation run, however, suggest that the modeling principles as well as their implementation are sound. The results show that EUHFORIA is successful at simulating the large-scale heliospheric dynamics including complex solar wind streams interacting with multiple eruptions, making it a useful tool not only for space weather prediction purposes but also for scientific studies of the heliospheric plasma environment.

The results also point towards ways in which the current baseline model can potentially be improved. The two modeled CMEs that were launched further away from the central meridian possibly suggest that the density of the CME front at the flanks is lower than provided by the cone model. More event studies in particular using multi-spacecraft events are required in order to test whether this is in fact the case. Further, the results suggest that the modeled CMEs are too narrow, and capture better the shock and adiabatic compressed plasma in the downstream than the wake of the eruption. Modeling the plasma properties in the wake of the CMEs accurately is important in particular when considering multiple eruptions as the dynamics of CMEs erupting later are influenced by the prevailing disturbed heliospheric condition.

Several efforts for improving and moving EUHFORIA beyond the baseline are currently on-going. A major next goal is to employ better CME models that take into account the internal magnetic field of the CMEs, e.g., using a spheromak-based model (Kataoka et al., 2009; Shiota & Kataoka, 2016). The key problem in employing magnetized CME models in routine space weather predictions is the lack of a practical method to estimate its intrinsic (i.e., immediately after the onset of the eruption) magnetic field properties as well as the evolution of CMEs in the corona wherein significant rotation and deflection in both longitude and latitude can occur (Kilpua et al., 2009; Isavnin et al., 2014). Tackling these issues remains a significant topical research problem. Another major path is to improve upon the accuracy of the semi-empirical solar wind model in order to better predict, e.g., arrival times and durations of high speed streams. One potential avenue is to relax the assumption of a steady state solar wind solution in favor of a time-evolving method (Hayashi, 2013, Hayashi 2012; Shiota et al., 2014; Hayashi et al., 2015; Merkin et al., 2016). In addition to more extensive validation efforts, these topics will be addressed in upcoming work.

Supplementary Material

Movie of Fig. 5: the radial speed Access here

Movie of Fig. 5: the number density Access here

Acknowledgements

These results were obtained in the framework of the projects GOA/2015-014 (KU Leuven), G.0A23.16N (FWO-Vlaanderen), C 90347 (ESA Prodex), University of Helsinki three-year grant project 490162 and the SolMAG project (4100103) funded by the European Research Council (ERC) in the framework of the Horizon 2020 Research and Innovation Programme. For the computations we used the infrastructure of the VSC Flemish Supercomputer Center, funded by the Hercules foundation and the Flemish Government department EWI. The editor thanks Pete Riley and an anonymous referee for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 Schematic of the steps involved in constructing the semi-empirical coronal model that provides the boundary conditions for the heliospheric MHD model.
In the text

	Fig. 2 The phases of the heliosphere simulation with typical durations. Each run starts with a solar wind relaxation, followed by a phase in which previously observed CMEs are introduced. The forecast starts at t = 0 and is typically run for 5 days.
In the text

	Fig. 3 Radial component of the magnetic field as given by the coronal magnetic field (PFSS) model at r = R⊙.
In the text

	Fig. 4 Open and closed regions as determined by the coronal magnetic field model. Blue (red) pixels correspond to open field lines that continue in to the solar wind with a magnetic field pointing toward (away from) the Sun. Gray pixels indicate regions of closed magnetic topology, with lighter gray indicating a positive field at the foot point.
In the text

	Fig. 5 Snapshot of the MHD simulation at 03:03 UT on June 21, 2015. Top row shows the radial speed, while the bottom row shows the scaled number density. The left panels depict the solution in the heliographic equatorial plane, while the right panels show the meridional plane that includes Earth.
In the text

	Fig. 6 Same as Figure 5 but at time 11:01 UT on June 22, 2015.
In the text

	Fig. 7 Same as Figure 5 but at time 11:01 UT on June 24, 2015.
In the text

	Fig. 8 Same as Figure 5 but at time 11:03 UT on June 27, 2015.
In the text

	Fig. 9 Radial speed at the position of Earth as a function of time. The blue curve shows the data from the simulation (saved at 10 min cadence), while the red curve shows OMNI 5 min data.
In the text

	Fig. 10 Number density at the position of Earth as a function of time. The blue curve shows the data from the simulation (saved at 10 min cadence), while the red curve shows OMNI 5 min data.
In the text